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Abstract. 

The hyperspherical harmonic (HH) method has been widely applied in recent times 
to the study of the bound states, using the Rayleigh-Ritz variational principle, and of 
low-energy scattering processes, using the Kohn variational principle, of A — 3 and 4 
nuclear systems. When the wave function of the system is expanded over a sufficiently 
large set of HH basis functions, containing or not correlation factors, quite accurate 
results can be obtained for the observables of interest. In this paper, the main aspects 
of the method are discussed together with its application to the A = 3 and 4 nuclear 
bound and zero-energy scattering states. Results for a variety of nucleon-nucleon (NN) 
and three-nucleon (3N) local or non-local interactions are reported. In particular, NN 
and 3N interactions derived in the framework of the chiral effective field theory and NN 
potentials from which the high momentum components have been removed, as recently 
presented in the literature, are considered for the first time within the context of the HH 
method. The purpose of this paper is two-fold. First, to present a complete description 
of the HH method for bound and scattering states, including also detailed formulas 
for the computation of the matrix elements of the NN and 3N interactions. Second, to 
report accurate results for bound and zero-energy scattering states obtained with the 
most commonly used interaction models. These results can be useful for comparison 
with those obtained by other techniques and are a significant test for different future 
approaches to such problems. 



1. Introduction 

The study of nuclear physics started just in the first few years after the introduction of 
quantum mechanics. From the very beginning, the non-relativistic Schrodinger equation 
was introduced, written as 

H^> = (T + V)m = . (1) 
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Here E and \1/ are the energy and the wave function of the considered nuclear system, 
and H is the Hamiltonian operator, given as the sum of the kinetic (T) and potential (V) 
energy operators. Correspondingly, two important and difficult problems immediately 
were raised: the first one was the determination of the nuclear interaction, i.e. the 
potential energy operator, the second was the solution of Eq. (1) to calculate the wave 
function \1>. These two aspects are reviewed in the following sections. 

1.1. The Nuclear Interaction 

Significant progress on the problem of determining the nuclear interaction was made 
after the fundamental paper of Yukawa [1] and the subsequent experimental discovery 
of the pion. However, the very first field theoretical approaches focusing on pion- 
exchanges [2, 3] were not successful, due to the non-renormalizability of meson field 
theory. As a consequence, semi-phenomenological approaches were then pursued. They 
were based on the assumption that the nuclear interaction is mainly a nucleon-nucleon 
(NN) interaction, whose long-range part can be described as due to a one-pion exchange 
(OPE). For the medium- and shorter-range parts, empirical forms were used, with the 
inclusion, in some cases, of a hard-core repulsion [4, 5]. 

The discovery of heavy mesons gave rise to other attempts to describe the NN 
interaction as due to a OPE, heavy meson and multi-meson exchanges, with the short- 
range part still treated in a phenomenological way. Moreover, a vast theoretical effort 
was made to calculate the two-pion exchange contribution producing the intermediate 
range attraction of the NN force. The Paris [6], the Nijmegen [7] and the Bonn [8] 
models are the most representative ones of such effort. 

In recent times, accurate calculations on few-nucleon systems have clearly indicated 
that, to perform a meaningful comparison with the experimental data, the input NN 
potential must reproduce, with a x 2 /datum close to unity, the deuteron binding energy 
and NN scattering data up to the pion threshold, namely E lab < 350 MeV, Ei ab being the 
NN laboratory kinetic energy [9, 10, 11]. Consequently, a new generation of interaction 
models, as the so-called Nijmegen [9], Argonne v 18 (AV18) [10], and charge-dependent 
Bonn (CDBonn) [11] potentials, has been derived. These models explicitly include 
charge symmetry breaking (CSB) terms in the nuclear interaction in order to reproduce 
equally well the np and the pp data. Furthermore, they can be local, as the first 
two models, or non-local, as the last one. However, it has been shown that none of 
these models reproduces the A > 2 dynamics in a satisfactory way. For example, they 
strongly underpredict the A > 2 binding energies. A possible cure is the inclusion of 
three-nucleon (3N) interaction terms in the potential energy operator. The very first 
model for a 3N interaction was proposed by Fujita and Miyazawa [12] and was again 
based on meson theory. In fact, the 3N interaction was derived from the exchange of two 
pions between three nucleons, with the intermediate excitation of a A resonance. This 
two-pion-exchange 3N interaction mechanism is the starting point of the 3N interaction 
models derived by the Urbana group [13, 14]. The Urbana models are written as the sum 
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of the two-pion-exchange 3N interaction plus a phenomenological repulsive term. The 
strengths of the two contributions are adjusted to reproduce the triton binding energy 
and the nuclear matter equilibrium density, in conjunction with one of the Argonne NN 
potentials. In particular, the so-called Urbana IX 3N interaction [14] is often used with 
the Argonne AV18 NN interaction. A more recent and sophisticated 3N interaction 
model [15], still derived by the Urbana group, contains two-pion-exchange terms due 
to pion-nucleon scattering in S- and P-waves, three-pion-exchange terms due to ring 
diagrams with one A resonance in the intermediate states, and again a phenomenological 
repulsive term. The model has five parameters which are fitted to the light nuclei mass 
spectrum. This more recent model has a rather complicated operatorial structure and 
is currently object of study. 

Another family of models for the 3N interaction is known as the Tucson- 
Melbourne [16] (TM) potential, which arises from an off-mass-shell model for the pion- 
nucleon scattering based upon current algebra and a dispersion-theoretical axial vector 
amplitude dominated by the A resonance. The model contains monopole form factors, 
whose cutoff is adjusted to reproduce the triton binding energy. More recently, the model 
has been revisited within a chiral symmetry approach [17], and it has been demonstrated 
that the contact term present in the TM model should be dropped. This new TM 
potential, known as TM', has been subsequently readjusted [18]. The final operatorial 
structure coincides with that one given in the 3N Brazil interaction already derived 
many years ago [19]. 

The main difficulty associated with the aforementioned approach is the following. 
The potential models considered fit the experimental NN scattering data up to Ei ab < 
350 MeV, and the A = 3 nuclei binding energies, but they treat each one in its peculiar 
way the higher momentum contributions. More in general, the problem is the connection 
between these NN plus 3N potentials and the fundamental theory of strong interactions, 
namely quantum chromodynamics (QCD). The approach of Weinberg [20, 21], when 
applied to low-energy processes, suggests how the problem can be overcome, at least 
partially. The starting point is to consider the most general Lagrangian that incorporates 
the assumed symmetry principles, in particular the (broken) chiral symmetry of QCD, 
constructed in terms of the pion and nucleon fields, which are the effective degrees of 
freedom at low energy, and their covariant derivatives. The old-fashioned perturbation 
theory can be applied to NN scattering, resulting in an infinite number of Feynman 
diagrams. Weinberg has shown [21] that a systematic expansion of the nuclear potential 
can be made in terms of (Q/A) v , Q being the pion momentum, A 1 GeV the chiral 
breaking scale and v a non-negative integer number. For each value of u, i.e. the 
order of the expansion, the number of corresponding Feynman diagrams is finite and 
can be calculated. Such an approach is model independent and is known as the chiral 
perturbation theory (%PT). Following the first approach of Weinberg, a number of NN 
potentials has been obtained by many authors [22, 23, 24, 25, 26, 27, 28] increasing 
the order of the expansion. For zy = 1,2,... the notation NLO, next-to leading order, 
N2LO, next-to-next, and so on, is commonly used. 
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It is important to notice that %PT makes specific predictions for many-nucleon 
forces, too. Three-body forces, for instance, appear for v > 2, i.e. firstly at N2LO. This 
explains why the contributions of the 3N interaction to the description of the A > 2 
dynamics are rather small if compared with the ones arising from the NN interaction. 

More recently a new class of NN interactions has been obtained. With the purpose 
of eliminating from the semi-phenomenological high precision potentials the strong high 
momentum parts, the Hilbert space is separated into low and high momentum regions 
and the renormalization group method [29] is used to integrate out the high-momentum 
components above a cutoff A. In this way a model independent potential Vi ow -k can 
be obtained [30, 31, 32, 33, 34]. An important issue of this procedure is what value 
should be used for A. As already discussed in Refs. [30, 35], A should be large enough to 
include all the relevant degrees of freedom. For instance, if A < 2 m w , m n being the pion 
mass, the two-pion exchange part of the bare NN interaction is integrated out. As a 
consequence, the model-dependence of triton binding energy is remarkably reduced [35]. 
On the other hand, such mo del- dependence becomes very strong for A > 2.0 — 2.5 fur 1 , 
where the differences between the various semi-phenomenological NN potentials in their 
short- distance structure start to be important. 

1.2. Calculations on Three- and Four-Nucleon Systems 

Different methods can be used to study few-body bound and scattering states. In nuclear 
physics many calculations have been made on systems with A = 3 and 4 by solving the 
Schrodinger equation using various NN local or non-local interactions with and without 
the inclusion of 3N forces. For A = 3, the Faddeev equations [36] (FE) technique has 
been widely applied both in coordinate and in momentum space. A detailed discussion 
of the method and the results for bound and scattering states can be found in Ref. [37]. 
The Faddeev- Yakubovsky (FY) equations [38] are the generalization of the Faddeev 
equations to the A = 4 case. They can be solved in configuration space [39] as well as 
in momentum space [40, 41, 42]. 

Alternative approaches to the FE and FY equations techniques make use of 
variational principles. Among them we can mention the Rayleigh-Ritz variational 
principle used for studying bound states and the Kohn variational principle [43] used 
in the case of scattering states. A few of the most important variational approaches for 
studying few-nucleon systems are briefly recalled here. 

The quantum Monte Carlo variational methods have been applied in two different 
ways. In the first one, known as the variational Monte Carlo, flexible functions 
containing a number of trial parameters are conveniently chosen. The mean value of the 
Hamiltonian is calculated by a Monte Carlo technique and the parameters are varied so 
as to obtain the minimum mean value [44]. The second approach is the so-called Green's 
function Monte Carlo method. It is based on a sophisticated algorithm for evaluating 
appropriate path integrals in order to obtain the properties of the bound states of the 
Hamiltonian [45, 46]. Further refinements of the algorithm have made it possible to 
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push the calculation up to A = 10 [47]. 

The coupled-rearrangement-channel Gaussian (CRCG) basis proposed by 
Kamimura [48] has been applied successfully to various systems, in particular to A = 3 
and 4 nuclei [49]. A further development of this approach is the so-called stochastic 
variational method (SVM) in which the Gaussian-type basis is optimized by means of 
a stochastic procedure [50, 51]. 

The no-core shell model (NCSM) is an approach which has been successfully applied 
both to s-shell and p-shell nuclei [52, 53, 54, 55, 56]. The calculation is done with a 
finite harmonic-oscillator basis. In the model space spanned by these states, an A- 
body operator is approximated to a two(three)-body effective interaction, by applying 
the Suzuki-Lee method [57], and therefore an effective Hamiltonian is constructed. By 
extending the basis, the NCSM calculation should converge to the exact solution. 

The hyperspherical harmonic (HH) method has a rather long history. In principle 
the HH basis functions can be used for expanding the wave function of a generic A- 
body system [58, 59] . The problem with such a basis is its quite large degeneracy which 
increases very rapidly with A and makes necessary the inclusion of a very high number 
of basis elements in the expansion. Exceptions to this general rule arise when the system 
interacts through very soft interactions or when the interaction has an hyperradial 
character as the harmonic oscillator potential. However, this is not the case of the 
strong state-dependent nuclear potential making the application of the HH technique to 
nuclear systems a difficult task. One possibility to tackle this problem is to introduce 
an effective interaction in the Hilbert space spanned by the finite basis considered [60], 
as has been done in the applications of the NCSM mentioned above. 

During the last few years, the authors of this paper started a strong collaboration 
in the application of the HH technique to describe bound states of light nuclear systems 
and to low-energy scattering processes involving few nucleons. In the first approaches, 
initiated by some of the authors, the HH functions were multiplied by suitable 
correlation factors in order to strongly improve the convergence of the basis [61, 62, 63]. 
Subsequently, thanks to the improvement in the computational facilities, it has been 
possible to consider very large sets of basis functions in order to accurately describe 
all the details of the wave function of the system [64, 65]. A significant test of the 
accuracy reached in the calculation of the wave function of the a-particle bound state 
using the HH expansion is given in Ref. [66]. In that work, the HH results for the 
binding energy and other quantities of interest are compared to those ones obtained 
by the methods quoted above. In the production of this benchmark, the Argonne v' 8 
(AV8') NN potential [45] has been used. Very recently, the HH technique has been 
further developed making it possible to use for the first time non-local potentials in the 
description of the A = 3, 4 systems [67, 68]. 

The HH method is a powerful tool for investigating scattering processes in the 
A = 3 [69, 70, 71] and 4 [72] systems. The scattering wave function is expanded in the 
HH basis taking into account explicitly its asymptotic part describing the relative motion 
of the incident nucleon and the target. The Kohn variational principle [43] is then used 
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to calculate the S-matrix of the reaction from which the observables of interest can be 
computed. To be noticed that in the application to scattering states the electromagnetic 
long range potential has to be taken into account properly. A discussion of this subject 
can be found in Refs. [73, 74]. Furthermore, different benchmarks have been done in 
which the rial and pal phases calculated using the HH expansion are compared to those 
obtained using the FE technique [75, 76]. As a further application, the scattering wave 
functions obtained from the HH expansion can be used to describe electroweak few- 
nucleon reactions. It should be remarked that at small energies, as those ones involved 
in astrophysical processes, the HH technique is the only one presently available that is 
capable of yielding precise predictions for a few astrophysical reactions [77, 78, 79, 80]. 

The present paper is organized in the following way: in Sec. 2 the Jacobi coordinate 
systems are introduced, together with the hyperspherical variables. In Sec. 3 the 
hyperspherical harmonic functions are defined and their properties are briefly reviewed. 
In Sec. 4, the practical implementation of the HH expansion for the A = 3 and 4 bound 
states is described. In particular, the two "standard" approaches known as HH and 
correlated-HH (CHH) expansions are discussed. In Sec. 5, the HH method devised for 
the study of low-energy scattering states is briefly reviewed. In Sec. 6 a large variety of 
results obtained applying the HH method to the A = 3 and 4 bound and scattering states 
is presented and compared with those obtained by other techniques and the available 
experimental data. Some concluding remarks, with an outlook on the future of this 
technique, are presented in Sec. 7. Finally, a few details of the calculations are given in 
the Appendices. 

2. Jacobi Coordinates and Hyperspherical Variables 

Let us consider a generic isolated system of A particles with spatial coordinates and 
masses wij, % — 1, . . . , A. In the center-of-mass reference frame, N = A — 1 vectors are 
sufficient to specify the spatial configuration of the system. They can be taken as linear 
combinations of the r^ and there is considerable freedom in choosing the corresponding 
coefficients. Let us introduce the Jacobi vectors Xj, i — 1, . . . , N, which, by definition, 
are such that the total kinetic energy operator can be written in the form 



where m is a reference mass, M = Yli=i m « * s the total mass of the system, and 






(3) 



where 




(4) 
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i.e. Xj is the coordinate of the center-of-mass of particles In this paper we will 

be interested in cases where all the particles have equal mass m and correspondingly 
Eq. (3) takes the form 



x JV-j+l 



2j 



J + l 



N . (5) 



For A — 3, the above definition of the Jacobi coordinates, once a permutation of (1, 2, 3) 
has been fixed, is unique. For A = 4 two different choices of Jacobi vectors exist, 
specifically 

set A set B 

X 1A — \ H r m 3 ) ' X l-B = r m ~ r fc , 

r r+r 4. r + r (6 



X 3A — r j — r i j X 3B — r j ~ r i 5 

where (i,j,k,m) is a generic permutation of (1,2,3,4). In the following, the wave 
function of the system will be expanded into a complete basis. The basis elements can 
be defined using set A or set B of the Jacobi coordinates. The completeness of the 
basis assures that the expansions using one set or the other are completely equivalent. 
In numerical applications the expansion is truncated and this equivalence does not hold 
anymore. Therefore it could be convenient to include basis elements defined both in 
set A and set B. As an important example, let us consider the ground state of the 
a-particle. The set A of coordinates in Eq. (6) can be more adequate for constructing 
contributions to the wave function corresponding to a {3 + 1} cluster structure, namely 
3 He + n or 3 H + p. Set B, instead, is more suitable for contributions coming from the 
{2 + 2} configuration, such as a d+d cluster structure. Therefore the use of an expansion 
basis constructed using both coordinate sets can speed up the convergence [49, 63]. 

For a given choice of the Jacobi vectors, the hyperspherical coordinates are given by 
the so-called hyperradius p, which results to be independent of the permutation order 
of the particles and is defined as 

P 2 = f> 2 = f 5> " r ^ = 2 I> " X) 2 , (7) 

i=l i<j i=l 

and by a set Qn of angular variables. In the Zernike and Brinkman [58, 81] 
representation, they are 2N polar angles x,i = (9 i: fa) of the Jacobi vectors Xj, 
i = 1, ... ,7V, and TV — 1 hyperspherical angles ipi, with < (pi < ir/2, given by the 
relation 

cosifi = —= 1 , i = 2,...,N, (8) 

yjx{ + . . . + Xf 

where Xi is the modulus of the Jacobi vector Xj. Therefore, we have 

fijv = {xi, ■ ■ ■ ,x N ,ip 2 , ■ ■ ■ ,<Pn} ■ (9) 
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3. Hyperspherical Harmonic Functions 

The spherical harmonic functions Y^ m (0, <p) are a very familiar tool in the study of 
quantum mechanical problems in a three-dimensional space. It is well known that 
T^Yi jTn ,(6, (ft) is a harmonic polynomial of order i in the Cartesian components x,y,z of 
the vector r, a harmonic polynomial being defined to be a homogeneous polynomial 
satisfying the three-dimensional Laplace equation. The spherical harmonic functions 
form an irreducible representation of the SO (3) group. The hyperspherical harmonic 
functions are the generalization to the case of D — 3N dimensional space. To this end, 
it is convenient to consider a homogeneous polynomial h[G] of degree G in the Cartesian 
coordinates of N Jacobi vectors and then introduce 

Y [G] = h [G] /p G , (10) 

where p is the hyperradius as given in Eq. (7). Since h\Q\ is a homogeneous polynomial, 
h[G]/p G does not depend anymore on p and therefore is a function only of the variables 
Qn specified in Eq. (9). Moreover, the 3iV-dimensional Laplace operator A can be 
written in the form 

where is the 3iV-dimensional generalized angular momentum operator depending 
only on the hyperangles Qjy. By definition Ah[c] = 0, so that 

Ah -f 92 i 3N ~ 1 9 i *Mh*y \ n a Y (D \ 

Ah[G] -\W ~~P~dp + ~^) p Y[g]{Qn) 

= (a 2 n (Q n ) + G(G + D- 2))p G " 2 F [G] (fi JV ) = , (12) 
and dividing by p G ~ 2 we get 

n) — . (13) 

A function satisfying the latter equation is an eigenfunction of the generalized angular 
momentum operator and is known as a hyperspherical harmonic function. The operator 
A 2 N is sometimes called the grand angular momentum operator, and G the grand angular 
momentum quantum number. [G] stands for a set of quantum numbers as it will be 
specified in the following. The grand angular momentum operator can be written in the 
form [58] 

d 2 r j Q 

A 2 (Vti) = t;— 2 + 3(i — 2)cotany?j + 2(cotan9?i — tany?/ 



dip 2 L J 

+ + A tf , (14) 

cos"' (fi sm ipi 

where — L 2 is the angular momentum operator associated with the i-th Jacobi vector. 
In particular 

Aim = L\ . (15) 
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The solutions of Eq. (13) can be constructed by following a recursive procedure [81]. In 
the case N — 1, namely D = 3, the eigenfunctions are the spherical harmonic functions 
Y timi (xi). 

For N = 2, the equation to be solved is 

A 2 (n 2 )Y [G] (n 2 ) = -g(g + a)y [g] (q 2 ) , n 2 = { Xl ,x 2 , V2 } , (16) 

with A 2 (£l 2 ) following from Eq. (14). Let us look for a solution of the previous equation 
of the form 

Y [G] (Q 2 ) = F(cos2 V 92)(cosv? 2 ) £2 (sin^ 2 )^r, imi (a; 1 )F £2m2 (x 2 ) , (17) 

where F is a function to be determined. In terms of the variable z = cos 2(p 2 , from 
Eq. (16) one gets 

(1- z 2 )F" + (a- [3z)F' + -fF = , (18) 

where 

a = e 2 -£i, P = ti + £ 2 + 3 , 

7 = \ [G(G + 4) - (£ 1 + i 2 ){i x +£ 2 + 4)] . (19) 

Eq. (18) is the one satisfied by the Jacobi polynomial Pn 1+l ^ 2,t2+l ^ 2 (z) provided that 
G = 2n + t\ + i 2 [82]. Therefore, a solution of Eq. (16) is given by 

Y [G] (n 2 ) = Af^(cos V2 y>(sm V2 )^Y eimi ( Xl )Y e2m2 (x 2 ) 

x P^ +1 / 2A+1 / 2 (cos2v9 2 ) , (20) 

where M^ ,v ' 2 is a normalization factor which will be specified later and v 2 = 2n + i x + 
£ 2 + 2. It is easy to verify that p G Y[ G ](Q 2 ) is a homogeneous polynomial of order G 
in the Cartesian components of X!,x 2 . In fact, for N = 2, one has x\ = psirnp 2 and 
x 2 = pcosy?2, and therefore 

p G Y [G] (Q 2 ) =Nt U2 p 2n x[ 1 Y limi {x l ) 4 2 Y e2m2 (x 2 ) 

x Pt +1/2/2+1/2 (cos2<p 2 ) . (21) 
Since p^ 1+1 / 2 ^ 2+1 / 2 j s a polynomial of degree n in the variable cos2y9 2 and 

p 2 cos 2y? 2 = 2x 2 2 - p 2 = x\ - x\ , (22) 

one has 

n 

p 2n p . 1+ l/2,£ 2+ l/2 (cos2v92) = p 2„ ^ Qm(cos 2( ^)™ 

m=0 

= ^a m (p 2 cos 2v, 2 )V ( ^ m) = ^a m (x 2 -x 2 r(a;i + x 2 r- , (23) 

m=0 m=0 

which is a homogeneous polynomial of degree 2n. Moreover the two terms in Eq. (21), 
x^Y timi (xi) and x^Y^JyX-i), are homogeneous polynomials of degree i\ and i 2 
respectively, so that in conclusion, p G Y^{VL 2 ) is a homogeneous polynomial of degree 
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G = 2n + l x + £ 2 . In this case the symbol [G] stands for the set of quantum numbers 
£i, mi, £ 2 , m 2 and n. 

From the expression of the grand angular momentum operator Af given in Eq. (14), 
it is possible to determine the HH functions for a given N in terms of those found for 
the case N — 1. Let us consider the function Y[c JV „ 1 ](fi7v-i), satisfying the equation 

A^Qiv-OY^.^Ar-i) = -G*_i(GW-i +D- 5)Y [Gtf _ l] (Q N . 1 ) , (24) 

namely the eigenfunction of A^_ x whose grand angular momentum quantum number is 
Gn-i, and look for the eigenfunction of of the form 

Y [G] (n N ) = (cos VN ) eN (sm VN ) G ^ 

x Y [GN _ l] (n N - 1 )Y eNmN (x N )F(cos2(p N ) . (25) 

By inserting this expression for Yjq in the corresponding eigenvalue equation and taking 
into account Eq. (24), the following solution is found 

F(cos2 VN ) = K n ' Un PZ- i/n+1/2 ( co ^n) , 



(26) 



with 



G =2n N + £ N + G N - 1 
1 . 



v N -i = G 



N-l 



+ 



3(N-1) 



(27) 



The complete expression of the HH function can therefore be cast in the form [58] 





' N 




Y [G] (n N ) = 




X 









JV 



J'=2 



where 



cos2<^ 



(28) 



(29) 



and the quantum numbers Gj and ia- are defined to be 



Gj = + 2n ») ' ni = ' G = Gn > ^ = G i + 4 - 1 • ( 3 °) 



i=i 



The symbol [G] distinguishes between different HH functions having the same G value. 
There are 3N — 1 quantum numbers which specify a HH function: the 2N quantum 
numbers £, m associated with the spherical harmonic functions, and the N — 1 quantum 
numbers n associated with the hyperspherical polynomials. In summary, [G] stands for 
the following set of quantum numbers 

[G] = {£1, ...,£ N , mi, ... , m N , n 2 , . . . , n N } . (31) 
The normalization factors M are chosen so as to verify the orthonormality condition 



dQ N (Y[ G ](Q N )^ Y[ G >](Q N ) = S[ G ],[g'] , 



(32) 
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(33) 



i=2 



is the surface element on the hypersphere of unit hyperradius. Their explicit expression 
is 

2u j Y(u j -n j )n j \ i 1 / 2 



M J 



Tin 



(34) 



Again, by using the definition of the hyperspherical angles in Eq. (8), it can be easily 
shown that the functions p G Y^(fl N ) are homogeneous polynomials. In fact, 



P G Y [G] (n N ) oc 



JV 



Ylx^Y t . tm .(xj) 



X 



' AT 



Yl P 2n > (sin <p j+1 ■■■sin (p N ) 2nj P nj (cos 2<pj) 

J =2 

The first factor is a homogeneous polynomial of degree t\ + h £n- Moreover, 



p cos2ipj = 



x 2 -(x 2 + --- + x 2 _ 1 ) 
(sinipj + i ■ ■ ■ sinip N ) 2 ' 



and 



Therefore 



(psmipj+i ■ ■ -sin^jv) 2 = x\ + V x 2 . 

(p sin tp j+1 • • • sin ^ N ) 2n] P nj (cos 2tpj) 

rij 

= p 2rij (sin(/9j + i • • •sin(^Ar) 2nj a m (cos 2(pj) r 



(35) 



(36) 



(37) 



m=0 



m=0 



(x 2 + --- + x 2 _ 1 )J [x\ 



+ --- + X 



0' 



(38) 



is a homogeneous polynomial of order 2rij. 

The HH functions have several important properties. A few of them will be reported 
here without proof. 



• The expansion of a plane wave in the 3N dimensional space is given by [59] 



where Q and Jl^ are the hyperspherical coordinates associated with the N Jacobi 
conjugate momenta q^, and J v (Qp) are Bessel functions of the first kind. The sum 
is taken over all the grand angular momentum quantum numbers G and over the 
corresponding degeneracy, namely the number of different HH functions for each 
G. 
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• By introducing the hypercoordinates p', Q' N associated with the N vectors x- and 
using the plane wave expansion given in Eq. (39), the following relation is found 

E^](^)%](Ov)^^ = f[5 3 (^-^ . (40) 

[G] ' i=l 

With the definition 

s D -\n N -a N ) = Y[ 5\x t - x>) n . d % . ^ 4 , (4i) 

we arrive at the relation 

^2 Y [G](^'n) y IG](^n) = S D -\n N - Q' N ) , (42) 

[G] 

which assures the completeness of the HH basis. As a consequence, every "regular" 
function /(^at) can be expanded in terms of the HH functions. In fact 

f(Q N ) = ( dtt' N 5 D -\Q N - n' N )f{n' N ) = J2aiG]YiG](tt N ) , (43) 

J IG] 

where 

a m = Jdn' N Y^(n' N )f(n' N ) . (44) 

• The 3iV-dimensional Fourier transform of a function /(xi, . . . , xjv) is given by 

/(qi, . . . , qjv) = ^ D/2 J d 3x i ■ ■ ■ d 3x N /(xi, . . . , x N ) 

x e i ^= 1 - Arqi ' Xi , (45) 

and the inverse transform is 

/(xi, . . . , xjv) = ^ D/2 J rf3{ li ' ' ' d3( lN /(qi, • • • , qw) 

x e"*^ 1 ^*'* 4 . (46) 

Note that 



A' 



-Lj / d^ ■ ■ ■ <*V e^- -(x,-x D = "Q 5 3 (Xi _ xD (4?) 

From the previous discussion it follows that the spatial wave function of an A-body 
system can be written in terms of the N Jacobi vectors and in turn expanded in 
the HH basis 

*(x 1 ,...,x JV )=^u [ G|(p)y[G](n JV ) , (48) 

[G] 

where the expansion coefficients are now functions of the hyperradius. The same 
expansion can be expressed as well in terms of the momentum space coordinates: 

9(q 1 ,...,q N )=^ m (Q)Y l a i (n q N ) ■ (49) 

[G] 
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The two expansions are equivalent provided that 

roo D-l 

9[G\(Q) = J d P (g p )£>/2-i j g+§-i(Qp) u IG](p) ■ (50) 

In this paper, HH functions with definite total angular momentum L, M are used. 
They are constructed (in configuration-space as an example) by the following coupling 
scheme, 

{A) H {G} , LM (n N ) = Yl (hm 1 £ 2 m 2 \L 2 M 2 )(L 2 M 2 e 3 m 3 \L 3 M 3 ) x 

mi,..., TUN 

• • • x (L N _ 1 M N _ 1 e N m N \LM)Y [G] (Q N ) , (51) 

where Y[g](Qn) is defined by Eq. (28), (£imi£ 2 m 2 \L 2 M 2 ) and so on are Clebsch-Gordan 
coefficients and Mj = Ylj=i i m r The symbol {G} of the function ^TC{g},lm(^n) now 
stands for the following set of quantum numbers 

{G} = {4, ■■■An, L 2 ,..., Ljv-i, n 2 , ■ ■ ■ , n N } . (52) 
The explicit expression of a HH function for A = 3 is 

= \ Y ^i)Y^2)] LM {2) Vi^{ V2 ) , (53) 

whereas for A = 4, 

(4) ^ lfe£ 3 L2n2n3}iL M(fi3) = [^0ri)^ 2 (£2)] i2 Yg 3 (£3) lm 

x (2) ^/ 2 (^2) (3) P^ 2+£l+£2/3 (^) • (54) 



4. HH and CHH Expansions 

In this section, a brief description of the HH expansion with or without the inclusion 
of correlation factors is presented. In the first subsection, we deal with the "standard" 
HH expansion for the systems of A = 3 and 4 nucleons. The transformation coefficients 
relating HH functions given in different sets of Jacobi coordinates are introduced in 
Sec. 4.2. It is well known that for realistic NN interactions the HH expansion is slowly 
convergent. An extensively used method to overcome such an unpleasant behavior is 
presented in Sec. 4.3, where the correlated HH expansion is introduced. The calculation 
of the correlation functions is discussed in Appendix C. In Sec. 4.4, some technical 
details of the practical applications of these expansions are reported. 



4-1. The HH Expansion 

Let us consider first the bound state of a trinucleon system with total angular momentum 
J, J z , and parity ir. The complete antisymmetrization of the state is guaranteed by 
writing the wave function as a sum of three Faddeev-like amplitudes 

* 3 = j2^ ( ?\4 ) ), (55) 
P =i 
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where 

*?> ='i~ *i , = ^ (r* - ^) , (56) 

are the Jacobi vectors corresponding to the p-th even permutation of the three particles 
(corresponding hereafter to the ordering i, j, k of the particles). The generic amplitude 
^(x^, x^) is, by construction, antisymmetric with respect to the exchange of particles 
i and j. Each ^(x^, x^) can be expanded in HH functions as 

N c N a M 

lKx?\ x«) = E E C (p) , (57) 

a=l n=N" 1=0 

where fi(p) is a complete set of functions for the variable p whose choice will be discussed 
in Sec. 4.4, and 

5™(p) = \ {3) n {eia e 2an} , La (n?) [S m 8 k ]A [T aa t k } TaTz . (58) 
Here is the HH function defined in Eq. (53) with 

(p) 

ee {x?,x®, , cos <p<?> = ^- , (59) 

p 

and Sk (tk) denotes the spin (isospin) state of particle k, whereas S aa (T aa ) the spin 
(isospin) of the i,j pair. Moreover, the total spin S a is coupled with L a to give the total 
angular momentum JJ Z . The set of quantum numbers 

a = {£i a , t2ai L a , S aa , S a , T aa , T a } , (60) 

is usually denoted as "channel" . The sum over a is truncated after the inclusion of the 
N c most important channels. Note that we allow for the possibility to have states of 
different total isospin T a (= 1/2 or 3/2) in our expansion. The quantum numbers of 
any channel must satisfy two requirements: in order to ensure the antisymmetry of the 
wave function, the amplitudes ip have to change sign under the exchange of the particles 
% and j, and therefore the number £ 2a + S aa + T aa must be odd. In addition, £ la + £ 2a 
must be an even or odd number in correspondence to the positive or negative parity of 
the considered state. The quantum numbers of the most important channels are given 
in Table 1. 

An important property of the basis elements is the following: due to the sum over 
the permutations in Eq. (55), some states 

C = t B n» (61) 
P =i 

are linearly dependent on others. Such states can be removed from the expansion basis 
by choosing appropriate values for the quantities N® in Eq. (57). As a matter of fact, 
the minimum to be considered depends on the set of channels and the HH functions 
included in the expansion of the wave function. For example, the states £>™ of the first 
two channels reported in Table 1 with n — 0, 1 are identical. This means that, if the first 
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Table 1. Quantum numbers and N® values for the first 14 channels considered in the 
expansion of the even parity and J = 1/2 wave function of a three- nucleon system. 
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two states of the channel a = 1 are included in the calculation, we must use iVi? = 2. A 
similar analysis has been made for all the channels and the minimum values N® needed 
to avoid linearly dependent elements have been reported in Table 1. The number of HH 
functions included in each channel is simply given by M a = N a — N® + 1. 

As a result of a number of investigations on the problem [64, 83, 84], it is well 
known that the components of ^(xj^, x^) with definite values of the angular momenta 
t\on d-2a give contributions rapidly decreasing when these values increase. As a matter 
of fact, the components with t la + £ 2 a > 6 give nearly negligible contributions in the 
case of all the available two-body NN interactions. On the contrary, very high values of 
the order n of the Jacobi polynomials must be considered in Eq. (57). 

We now generalize the HH approach to the four-body system. The wave function 
of a A = 4 system having total angular momentum J, J z , and parity it can be written 
as 

*4 = J2 [M^l x?j, x?j) + M^l xg, x?i)] , (62) 
P =i 

where the sum is taken over the twelve even permutations p of the particles. In the 
remainder, we suppose that the permutation p corresponds to the order i, j, k, m of 
the particles. The vectors x^, x^, x^ and x^, x^, x^ are the two possible sets of 
Jacobi vectors defined in Eq. (6). Each amplitude ip is then expanded in terms of the 
HH functions defined in Eq.(54): 

«x&,x£>,xW) = £ E E C *°.W MP) BZmsaip) , (63) 

a=l n2,«3=0 1=0 
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where X = A or B and 

nHH / \ 



^B,n 2 n 3 a(p) 



x [[Ms™ 5 * 

{ (4)/ ^laW3«£2a 



},L a (^' P) )x 



ri2n 3 },L 



^6o 



T a T z 



X 



- Sn 



(n?*)) x 



, (64) 



(65) 



Here, f^'^ denotes the hyperangular variables constructed with the set X of the Jacobi 
vectors relative to the permutation p of the particles. Moreover, (tj) denotes the spin 
(isospin) state of particle % and 

Oi = {^lo, ^2a> ^3a, L201, L a , S aa , Sba, S a , T aa , Tb a , T Q } , (66) 

specifies the generic four-body channel. In Eqs. (64) and (65) L a and S a are coupled to 
give J, J z . In order to ensure the antisymmetry of the wave function, the amplitudes 
ipA and ipB have to change sign under the exchange of the particles % and j. Therefore, 



the number 



• 3q 



+ s a 



T nn must be odd. In addition, 



'la 



-2a 



+ 



- 3a 



must be an even 



or odd number in correspondence to the positive or negative parity of the considered 
state. A few of the channels taken into account in the calculations are given in Table 2. 

As in the three-body case, the contribution of channels with l\ a + £ 2a + ha > 6 can 
be disregarded. However, the sum over n 2 and n 3 usually includes a very large number 
of terms. It has been found convenient to limit the sum over n 2 and n 3 in Eq. (63) to 
include HH functions with < n 2 + n 3 < N a , N a being a non- negative integer. For 
a given N a , the number M a of functions included in the expansion of the channel a is 
given by 

M a = (N a + l)(N a + 2)/2 . (67) 
After antisymmetrization of the state, some of the functions 



12 



p=l 



(P) 



(68) 



are linearly dependent on others and they must be removed from the expansion. As an 
example, it can be verified that £>™„ 3a=1 = £>™n 3 a=2 for all n 3 values (the channels 
a = 1,2 are specified in Table 2). In this respect, it is useful to know how many 
linearly independent states B exist for a given set of values of G, L, S and T, with 
G = £i a + £ 2 a + ^3a + 2n 2 + 2n 3 . Such a number has been calculated in Ref. [85], where 
it has been shown that it is significantly less than the degeneracy of the basis. As a 
result, a large number of states has to be removed from the expansion and, after that, 
the number M a results to be appreciably lower than that given in Eq. (67). 



A High- Precision Variational Approach 



17 



Table 2. Quantum numbers for the first 23 channels considered in the expansion of 
the even-parity and J — wave function of a four-nucleon system for the set A of 
Jacobi coordinates. 
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4-2. The Transformation Coefficients 

As discussed previously, in the case of A identical fermions, antisymmetrical states can 
be obtained by appropriate sums over the particle permutations. The calculation of 
the matrix elements of the potential between such states presents noticeable difficulties. 
In general, the integrals to be evaluated involve the potential and two HH functions 
constructed in terms of Jacobi vectors corresponding to different permutations of the 
particles. This task is simplified by the knowledge of the following transformation 
coefficients (TC), 

( %^(fift = E {A) -H {g% lm^n) , (69) 

{G'} 

where the sum is over all the quantum numbers {G 1 } specified in Eq. (52), with the 
condition that the grand angular momentum G defined in Eq. (30) is conserved, G' = G. 
As the number of functions with a given G value is finite, also the sum in Eq. (69) is over 
a finite number of terms. In Eq. (69), Q$ (fijv) specifies the hyperangular variables 
constructed with the permutation p of the particles. To simplify the notation for p = 1, 
in this section the subscript "1" will be omitted. 
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For the three-body system, these coefficients can be easily calculated using 
recurrence relations [86]. However, the TC can also be calculated by taking into account 
the orthonormality of the HH basis, namely 

a {G}%'} = J dQ » [ {A) ^{G'},lm(^ • (70) 

For A = 3, the above integrals reduce to bi-dimensional integrals of polynomial functions 
which can be calculated exactly using Gauss quadrature as discussed in Appendix A. 
For A = 4, Eq. (70) results in a 5-dimensional integration. In that case it is more useful 
to use a recurrence formula, derived in Ref. [85], and briefly outlined in Appendix B. 

The spin-isospin states appearing in Eqs. (58), (64) and (65) can be transformed 
in terms of spin-isospin states constructed using the reference order 1,2,3,... of the 
particles. The resulting states have the same total spin S a and total isospin T a as 
the initial states. Therefore, the complete functions B^ a {p) and -B™ 2n3Q (p) can be 
transformed to the reference permutation p = 1 and a given choice of Jacobi vectors. 
For example, for A = 4, 

£™ 2 na» = E 42T(p) B Z 2 n W (P = 1) , Ul = K, n> 3 , a'} ,(71) 



'2 "3 



where the sum is restricted to the indices n' 2 and n' 3 and channels a' such that 
G = £ la + £ 2a + ha + 2n 2 + 2n 3 = t la > + £ 2a > + ha' + 2n 2 + 2n' 3 , L a = L a >, S a = S a > 
and T a = T a i. Hence, summing over all the even permutations, it is possible to define a 
completely antisymmetric basis state given in the reference system p — 1: 

12 

ErHH /„\ — r}GL a S a T a ,Jir \ A AGL a S a T a ,Jn rjHH /„ t\ /70\ 

a X,ninzcXV) = &X,ii ~ 7 Ax.u.a' ^ 'A,n' 2 n' 3 a> \P ~ 1 ) > \ U ) 

p=l n' 

with jj, = {112, ri3, £i a , £ 2 a, ha, L 2a , S aa , Sb a , T aa , Tf, a } and 

12 

AGL a S a T a ,Jn _ \ y X,n2n- A a i \ (•70\ 
A X,/j,,ll> - / , a nLn' 3 a> VP)' \' 6 ) 

p=l 

The same analysis holds for the A = 3 system, where now /i = {n 2 , £i a , £ 2a , S aa , T aa } 
and p! = {n' 2 , a'}. Finally, the wave function of an A = 3, 4 system can be expanded as 

*r= E E^ifC/ T,j7r /Kp), (74) 

GLSTX n,l 

where the coefficients C^^ T are linear parameters (for A = 3, the sum over X is 
restricted to the single choice of Jacobi vectors specified in Eq. (56)). The above 
expansion can also be given in momentum space: 

E Y, C ^7®Z ST ' J * 9oAQ) , (75) 

GLSTX n,l 

where B ( ^ L ^ T ' Jn is now the hyperspherical-spin-isospin amplitude in momentum space 
and, according to Eq. (50), 

d P (g p )3JV/2-l J G+3N/2~l(Qp) fl(p) ■ (76) 
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Sometimes it will be convenient to consider HH functions having a specific value of 
j, the total angular momentum of particles 1,2. To this aim, it is useful to write the 
basis in the jj-coupling scheme 

gGLST,Jn _ ^GLST,Jn ^GTJn /yrA 

v 

The basis functions S^ TJ?r are HH functions in which the angular-spin part for A = 3 
is coupled as 

(78) 

JJz 



Ye 2 (x 2 )S a ) Y^i)^ 



and, for A = 4, as 

\\(Y e3 (x 3 )S a ) (Y i2 (x 2 )s 3 ) 1 (Y^sA } . (79) 

The coefficients Ax L ^J Jn are related to A^ L ^, j7T via Wigner coefficients. Now, the 
integer index v labels all possible choices of 

iy = {n 2 ,£2,S a ,j a ,e 1 ,j b ,T a } , forA = 3, (80) 
u= {n 3 ,£ 3 ,S a ,j 3 ,n 2 ,£ 2 ,j2,J2,h,ji,T a ,T b } , forA = 4, (81) 

compatible with the given values of G, T, J and n. 
4.3. The PEE and CEE Expansion 

The rate of convergence of the HH expansion results to be rather slow when the 
particle interaction contains large repulsion at small distances. For example, in the 
calculation of the trinucleon binding energy with the AV18 NN potential, HH functions 
with G up to 180 were considered [61]. The resulting expansion included a large 
number of terms (approximately 600). Only in that way, the calculated binding energy 
was found to be in agreement with the corresponding estimates obtained by other 
techniques [87, 88, 89, 90, 91]. The calculation rapidly becomes more and more involved 
when larger systems are considered. For the a-particle, for example, approximately 4000 
states have to be included in Eq. (62). 

In order to overcome this difficulty, in Ref. [92] the radial dependence of each 
amplitude has been modified by the inclusion of a suitably chosen correlation factor 
(correlated-HH or, briefly, CHH expansion). The role of the correlation factors is 
to speed up the convergence of the expansion by improving the description of the 
system when two particles are close to each other. In such configurations, there are 
large cancellations between the contributions from kinetic and potential energy terms 
and therefore the wave function must be very precisely constructed. Thus, it can be 
convenient to include from the very beginning in the wave function suitable terms for 
describing those configurations. In this way, the number of basis functions necessary to 
get convergence is strongly reduced. 
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For the three-body system, the functions B™(p) m Eq. (57) are replaced by the 
following ones, 

B^(P) = Fj$ B™(p) , (82) 

(3) 

where Fa P is a correlation factor which, in general, will depend on the channel a and 
the permutation p of the particles. In order to maintain the antisymmetry of the wave 

(3) 

function, Fa P is taken to be symmetric under the exchange of particles % and j. A choice 
of Fa P \ which has been extensively investigated, is 

F S = fa(rij)g a (r jk )g a (r ik ) , (83) 

where f a and g a are unidimensional functions determined as discussed in Appendix C. 
For a four-body system, the functions 

B^n 2n Up) = F S bZ.uUp) (84) 

are used to replace -B™ 2 n 3 a(p) in Eq. (63). In this case, the correlation factors have 
been taken of the form [63] 

F S = f a (nj)ga(r ik )g a (r jk )g a (r im )g a (r jm )h a (r km ) , (85) 

where f a , g a and h a are unidimensional functions of the interparticle distances. The 
choice of these functions is also discussed in Appendix C. 

The correlation factors so far introduced are of the Jastrow type, namely they are 
given by products over all the particle pairs. A simplified choice consists in correlating 
only one pair. The expansion basis obtained is called the pair-correlated hyperspherical 
harmonic (PHH) expansion, which, for a three-body system, is given by [62] 

B™{P) = /«(r y )S™(p) • (86) 

The use of this (simplified) correlation factor has been shown to provide very accurate 
results. 

It should be noticed that, with the inclusion of the correlation factors, the 
linear dependence between the states discussed in the previous subsections has to be 
reconsidered. In fact, for the A = 3 case, the minimum values N® of the index n entering 
in Eq. (57) can be taken to be iV° = for all the considered channels. Also for A = 4, 
the linear dependent states reduce to a small number (the actual number depends on 
the particular choice of the correlation factors). 

As for the HH expansion, it is possible to define a completely antisymmetric 
correlated state for A = 3 

3 

3L = £ S -(P)' (87) 
p=l 

and for A = 4 

12 

Bx,n 2 n 3 a = B X,n 2 n 3a (p) » (88) 

p=l 
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where A = CHH or PHH. Accordingly it is possible to expand the A = 3, 4 wave function 
in such a basis 

Ina 

E Cx„ 2 n 3 ^(p)£x,n 2 n 3Q . (89) 

A difference with the HH expansion is that the basis states are not labelled with the 
grand angular momentum G, due to the presence of the correlation factors. Furthermore, 
the presence of correlation factors complicates the transformation of the states given in 
permutation p to the reference system p — 1. As a consequence, multidimensional 
integrals have to be calculated to obtain the matrix elements when using the correlated 
basis states. This can be done very efficiently for A = 3 whereas, for A = 4, a stochastic 
integration method should be used. 



4-4- Details of the Calculation 

In general, the A = 3 and 4 bound state wave functions can be cast in the form 

W>=Z>I*C> , (90) 

c 

where |^) is one of the complete set of antisymmetric states defined in the preceding 
sections, and ( is an index denoting all the quantum numbers necessary to completely 
specify the basis elements. The coefficients of the expansion can be calculated using the 
Rayleigh-Ritz variational principle, which states that 

(6 C * J /\H-E\* J /) =0 , (91) 

where 5 c ^j^ indicates the variation of the wave function for arbitrary infinitesimal 
changes of the linear coefficients c^. The problem of determining and the energy E 
is then reduced to a generalized eigenvalue problem, 

E (*c I H ~ E I *C ) C C = ■ ( 92 ) 

C 

The main difficulty of the method is to compute the matrix elements of the Hamiltonian 
H with respect to the basis states |\E^). Usually H is given as a sum of terms (kinetic 
energy, two-body potential, etc.). The calculation of the matrix elements of some parts 
of H can be more conveniently performed in coordinate space, while for other parts it 
could be easier to work in momentum space. 

In the case of the HH expansion, the basis states are given in coordinate 
(momentum) space in Eq. (74) (Eq. (75)). The criteria to choose hyperradial functions 
fi(p) have been the following: (i) fi(p) — > for p — > oo; (ii) fi(p) should constitute an 
orthonormal basis; (iii) fi(p) should be easy to handle, when the Fourier transform of 
Eq. (50) is performed. A possible choice used here is given by 



flip) = 7 D/ %/ a + ^_ 1) , ^r X) (7P) e-*' , (93) 
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where l\ D 1 ' ) (7p) are Laguerre polynomials [82]. Here, there is only one non-linear 
parameter, 7, to be variationally optimized. In particular, 7 can be chosen in the 
interval 2.5-4.5 fm _1 for the AV18 and CDBonn potentials and 4-8 fm _1 for the N3LO- 
like potential, for both A = 3 and 4. The corresponding functions ga,i{Q) can be written 

as 

x 7: ^nrr P lfn/2 D '\u) , (94) 

(1 — u 2 ) 4 2 

where w = — ^=== = , P™ is an associated Legendre function and b l k is given by 

4 = nr('7^l- (95) 

so thatt< D - 1 '(i) = ELoH^ 

In such the normalization (iV) and kinetic energy (T) operator matrix 

elements can be easily computed, both in coordinate or in momentum space. Explicitly 

they are 

N%&%v = Effi J ')^% , (96) 
n" 

rpGLSTJn _ V^/ /.GLSTJttn* AGLSTJnrp / Q7 % 
1 X t d,X'ti'V — / J\ A X,nn" I A X>,n'n" ; l y 'J 



with 

= d Q Q 3A ~ 2 9g,v (Q) 9g,i (Q) , (98) 

where use has been made of the fact that the HH functions are eigenfunctions of the 
operator Ax(Q N ) defined in Eq. (13) corresponding to the eigenvalues G{G + 3A — 4), 
and that the functions fi(p) form an orthonormal set with respect to the weight 
p 3A ~ 4 . The calculation of the two-body potential energy matrix elements is more 
conveniently performed using the hyperspherical-angular-spin-isospin basis elements 
E^ TJn , which have a well defined angular momentum j between particles 1,2, and 
are totally antisymmetric. Therefore, the following relation holds 

< SjT^lSf * >= < 3gT'>(l, 2)|Sf * > . (99) 

The potential u(l,2) acts on the particle pair 1,2, accordingly to the quantum number 
j and all other quantum numbers in p, with the exception of S aa , T aa and £i jCe (£i, a /2,a) 
for A = 3 (A = 4), which are conserved. This reduces considerably the total number 
of matrix elements needed in the calculation. Moreover, using the TC introduced in 
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Sec. 4.2, the above equation reduces to a two-dimensional integral in one hyperangle 
(cos 03 for A = 4 or cos for A = 3) and p. 

In the case of non-local potentials, one has 

v(l,2) = V(x' N ,x N ) (100) 

in coordinate space, and 

v(l,2) = V(k> N ,k N ) (101) 

in momentum space, where the spin-isospin-dependence is understood. The integrals 
of Eq. (99) then are three-dimensional, and can be calculated using Tchebyshev 
and Laguerre weights and points for the hyperangular and the hyperradial variables, 
respectively [82]. A sufficiently dense grid can be used to obtain relative errors < 10~ 6 
for these integrals. More details about the procedure to compute matrix elements of the 
NN and 3N interactions are given in Appendix D. 

In the case of the CHH or PHH expansion, the matrix elements of the norm, 
kinetic energy and potential energy are all calculated numerically. These expansions 
have been used to treat the case of local potentials. The corresponding integrals reduce 
to three dimensions for A = 3 and to six dimensions for A = 4. In the first case a 
Gauss, Tchebyshev, Laguerre integration has been used for the variables pu — x\ • £2, 
z = cos 20, and p, whereas a quasirandom technique has been used for A = 4. 



5. The HH Technique for Scattering States 

We consider in this section the application of the HH expansion to a scattering problem. 

In particular, we focus our attention to elastic processes of the type N + Y — > N + Y, 

where N is a nucleon and Y is a bound nuclear system (A Y + 1 = A = 3, 4), in the low 

energy region, where the nucleus Y cannot be broken. 

The wave function ty^y* describing a N — Y scattering state with incoming orbital 

->-*-* 

angular momentum L and channel spin S (S = \ + Sy), parity 7r = (— ) L , and total 
angular momentum J, J z , can be written as 

m L N s ™* = m^ JJ * + m L A SJJ * , (102) 

T Q J J 

where ^ c z describes the system in the region where the particles are close to each 
other and their mutual interactions are strong, while fy^ SJJz describes the relative 
motion between the nucleon iV and the nucleus Y in the asymptotic region, where 
the N — Y nuclear interaction is negligible. The function -qjLSJj z ^ w \ 1 [ c \ 1 ^ Q van ish 
in the limit of large intercluster separations, can be expanded either in the HH or CHH 
basis as it has been done in the case of bound states (see Eqs. (55) and (62) for the 
HH and Eq. (89) for the CHH expansion). Therefore, applying Eq. (90), the function 
t^lsjJz can k e cas t ec [ j n the form 

i^c 5JJz ) = E c d^) > ( 103 ) 

C 

where |^) is defined in Sec. 4. 
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The function \fj^ SJJ z j s the appropriate asymptotic solution of the relative N — Y 
Schrodinger equation, and can be written as a linear combination of the following 
functions 

J2 \\Xn ® Ms ® Y L (y p j\ K x L (y p ) , (104) 



LSJJz /TV 



V 



where C is a normalization factor (see below). Here the sum over p has to be done over 
the even N p permutations of the A nucleons necessary to antisymmetrize the functions 
^lsjj z i (Xn ) is the Y (JV) wave function, y p is the distance between N and the center 
of mass of Y, Yi(y p ) is the standard spherical harmonic function, and the functions 
lZ\(y p ) are the regular (A = R) and irregular (A = /) radial solutions of the relative 
two-body N — Y Schrodinger equation without the nuclear interaction. These regular 
and irregular functions, denoted as J r L(y P ) and Ql^Vp) respectively, have the form 



(2L + l)\\q^C L ( V ) 



g L (y p ) = (2L + m L+l cMf R {y P )^f^- , (105) 

where q is the modulus of the N — Y relative momentum (related to the total kinetic 

2 

energy in the center of mass (cm.) system by T c m . = |-, // being the N — Y reduced 
mass), i] = 2fie 2 /q and £ = qy p are the usual Coulomb parameters, and the regular 
(irregular) Coulomb function Fl(t),£) (G\l(?7,£)) and the factor Cl(t)) is defined in the 
standard way [93]. The factors (2L + l)\\q L CL{rj) have been introduced so that T and 
Q have a well defined limit for q — > 0. The function fn{y p ) = [1 — exp(— by p )] 2L+1 
has been introduced to regularize Gl at small values of y p . The trial parameter b is 
determined by requiring that /i?(y) — > 1 for large values of y p , thus not modifying the 
asymptotic behaviour of the scattering wave function. A value of b ~ -Ry 1 , Ry being 
the dimension of nucleus Y, has been found appropriate. The non-Coulomb case of 
Eq. (105) is obtained in the limit e 2 — > 0. In this case F L (£)/£ and reduce to 

the regular and irregular Riccati-Bessel functions and the factor (2L + l)\\C L (r}) — > 1 
for 77 — > 0. With the above definitions, \l>^ 5JJz can be written in the form 

K s JJz = E [^^n^'jj. + Ks^AM's'jj] , (106) 

L'S" 

where the parameters ^^ 5 L , 5 ,(g) give the relative weight between the regular and the 
irregular components of the wave function. They are closely related to the reactance 
matrix (/C-matrix) elements which can be written as 

K-Uusi = (2L + 1)!!(2L' + mq L+L ' +1 C L (r))C L ,( V )TZ J LStL/sl . (107) 

By definition of the /C-matrix, its eigenvalues are tan 5lsj, 8lsj being the phase shifts. 
The sum over V and S' in Eq. (106) is over all values compatible with a given J and 
parity n. In particular, the sum over V is limited to include either even or odd values 
since (— 1) L ' = it. 
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The matrix elements TZ~[ S L'S'iv) an< ^ the l mear coefficients occurring in the 
expansion of \fjLSJJ z n are determined applying the Kohn variational principle [43] , which 
states that the functional 

[Ks,uAq)\ = Ks,uAq) - (tf£f/ J - \H - E\ K S -y) (108) 

has to be stationary with respect to variations of the trial parameters in ^^f-Y*- H ere E 
is the total energy of the system, and the normalization of the asymptotic states (factor 
C in Eq. (104) ) has been fixed by the condition: 

(n* SJJz \H - E\n[ SJJz ) - (nisjjjH - e\^ sjj j = 1 . (109) 

Using Eqs. (103) and (106), the variation of the diagonal functionals of Eq. (108) with 
respect to the linear parameters leads to the following system of linear equations: 

J2c C <V<\H - E\^ c >= -D x LSJJz (C) ■ (HO) 
C 

Two different inhomogeneous terms D x corresponding to A = R, I are introduced and 
are defined as 

DIsjjSO =< *t\H - E\n x LSJJz > . (Ill) 

The matrix elements T^'ls l's' are obtained varying the diagonal functionals of Eq. (108) 
with respect to them. This leads to the following set of algebraic equations 

^ T^LS,L"S"^L'S',L"S" = YLS,L>S' , (112) 
L"S" 

with the coefficients X and Y defined as 

X L S,L'S' =< ^LSJJ Z + ^c SJJzJ \ H ~ E \^L'S'JJ Z > > 

Y LS ,us> = -< n* SJJz + V L C SJJ - R \H - E\n[ ISIJJz > , (113) 
where \$/LSJJ z ,\ - g ^ e solution of the set of Eq. (110) with the corresponding 
inhomogeneous term. A second order estimate of r R-f jSL , s , is given by the quantities 
\n J LS us ,}, obtained by substituting in Eq. (108) the first order results. 

In the particular case of q = (zero-energy scattering), the scattering can occur 
only in the channel L = and the observables of interest are the scattering lengths. 
Within the present approach, they are defined as 

(2J+1 W = -lim^ J (114) 
and the corresponding asymptotic states are then given by 



xll 0JJJ z 
* A 



(115) 
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6. Results 

In this section results are presented for bound and zero-energy scattering states obtained 
within the HH approach, with and without the inclusion of correlation factors, using 
several interaction models. The section is divided into the following parts: in Sec. 6.1, 
the triton and 4 He bound states are studied with 5 different central potentials, which 
have been used by several groups to produce benchmark calculations. For some of 
them, different versions appear in the literature causing problems at the moment of 
comparisons. For this reason, the potential parameters are specified and the HH results 
are compared to those ones obtained with other techniques. In Sec. 6.2 the A = 3 and 
4 bound states are studied using the most recent models for the nuclear interaction, 
consisting of two- and three-nucleon potentials. Their predictions for binding energies 
and other ground state properties are presented. When possible, comparisons with other 
methods are also given. Finally, in Sec. 6.3, the results for the nd, pd, n 3 H, and p 3 He 
scattering lengths are presented and discussed. 

6.1. Bound States with Central Potentials 

Let us firstly consider the A = 3, 4 bound states using 5 central two-nucleon potential 
models, i.e., the Volkov [94], Afnan-Tang S3 (ATS3) [95], Minnesota [96], and Malfliet- 
Tjon version V (MT-V) and \/\\\ (MT-I/III) [97]. These potentials have been used by 
several groups to describe bound states of light nuclei. We review here briefly their main 
characteristics. Each potential V(r) can be written as 



where the function /(//j,r) is either exp(— /i^r 2 ) for the Gaussian-type potentials or 
exp(— /ijr)/r for the Yukawa-type potentials. The operators Vi act on the spin-isospin 
degrees of freedom and are written as 



where P r , P a and P T are the space-, spin- and isospin-exchange operators. The values 
of parameters V^, //j, Wi, Mi, B t and Hi of the potentials used in the present work 
are listed in Table 3. Few remarks are here in order: (i) the Volkov and MT-V are 
spin-independent models, while the other 3 potentials are spin dependent; (ii) the MT- 
I/III version acts only on s-waves; (iii) it is customary to include the point-Coulomb 
interaction (e 2 = 1.44 MeV fm) with the Minnesota potential; (iv) for all the considered 
potentials, the total orbital angular momentum is a good quantum number and therefore 
we have included in the wave functions only the channels with L = 0. Finally, note that 
the Volkov, ATS3, Minnesota and MT-V potentials are equal to those used in Ref. [98], 
while the version of the MT-I/III here adopted is the same as the one reported in Table I 
of Ref. [99]. 

The triton and 4 He binding energies have been calculated for these 5 central 
interaction models using the HH or the CHH approaches. The results are compared 




(116) 



Vi = ViX (Wi + M t P r + B t P a - H t P T ) , 



(117) 
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Table 3. List of the parameters of the central NN potentials used in this paper 
introduced in Eqs. (116) and (117). The potential strengths V, are in units of MeV 
for Gaussian- (G) and MeV fm for Yukawa-type (Y) potentials, respectively. The 
parameters Wi, Mi, Bi and Hi are dimensionless and the ranges in are in units of 
fm~ 2 for G or fm -1 for Y potentials, respectively. The Majorana mixture parameter 
M of the Volkov potential and the parameter u in the Minnesota potential are and 
1, respectively. 



Potential 


Type 


i 


Vi 


Mi 




M t 




Hi 


MT-V 


Y 


1 


1458.047 


3.11 


1.0 


0.0 


0.0 


0.0 


[97] 




2 


-578.089 


1.55 


1.0 


0.0 


0.0 


0.0 


MT-I/III 


Y 


1 


1438.72 


3.11 


1.0 


0.0 


0.0 


0.0 


[97] ' 




2 


-570.4255 


1.55 


1.0 


0.0 


0.0 


0.0 






3 


-56.4585 


1.55 


0.0 


0.0 


1.0 


0.0 


Volkov 


G 


1 


144.86 


0.82~ 2 


1.0 — M 


M 


0.0 


0.0 


[94] 




2 


-83.34 


1.60- 2 


1.0 — M 


M 


0.0 


0.0 


ATS3 


G 


1 


1000.0 


3.0 


1.0 


0.0 


0.0 


0.0 


[95] 




2 


-326.7 


1.05 


0.5 


0.0 


0.5 


0.0 






3 


-166.0 


0.80 


0.5 


0.0 


-0.5 


0.0 






4 


-43.0 


0.60 


0.5 


0.0 


0.5 


0.0 






5 


-23.0 


0.40 


0.5 


0.0 


-0.5 


0.0 


Minnesota 


G 


1 


200.0 


1.487 


0.5m 


1.0 - 0.5u 


0.0 


0.0 


[96] 




2 


-178.0 


0.639 


0.25u 


0.5-0.25u 


0.25u 


0.5 - 0.25u 






3 


-91.85 


0.465 


0.25m 


0.5 - 0.25u 


-0.25u 


-0.5 + 0.25u 



to those ones obtained by other techniques in Table 4. In all of these calculation, 
h 2 /m = 41.47 MeV fm 2 . Some of the techniques considered in the comparison have 
already been described in Sec. 1, and they are the FY equations method [40], the 
SVM [50, 51] and the CRCG expansion technique [49]. Furthermore, the results of 
the effective-interaction hyperspherical harmonics (EIHH) method of Ref. [100] and the 
ATMS (Amalgamation of Two-body correlation into the Multiple Scattering process) 
method of Ref. [101] have been reported. 

From inspection of the table, we can conclude that: (i) for the Volkov potential 
with Majorana parameter M — 0, our results agree very well with the estimates of the 
SVM [98] method. Note that the Volkov potential, given as a sum of Gaussians, has a 
very soft core and therefore the induced two-body correlations in the ground state wave 
function are weaker than in the other cases. In fact, we have found that the convergence 
of the HH expansion in this case is quite fast. Since only the inclusion of HH states 
with fairly low values of the grand angular momentum quantum number is sufficient to 
obtain convergence, a successful HH calculation for this potential was already possible 
more than 30 years ago for A = 3 [102] and more than 20 years ago for A = 4 [103]. (ii) 
Both for the ATS3 and Minnesota (with exchange parameter u—1) potential models, 
there is a good agreement between the different theoretical estimates. Note that both 
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Table 4. The triton and 4 Hc binding energies in MeV, calculated for various central 
interaction models, are compared with the results obtained by other methods. 

Potential Method B( 3 H) B( 4 He) 



Volkov 


HH 


8.465 


30.420 




SVM 


8.46 


30.424 


Al DO 


rirl 


O. / Do 


ol .Olo 




PIT IT 

V III! 


& 7^8 


Ol.Ol 




G.VM 
o V 1V1 


o. / Do 


ol .010 


Minnesota 


HH 


8.3858 


29.947 




V I 1 I 1 


O.OOOO 


zy.yo 




O V 1V1 


o.ooU 


on qq7 




T71TTTTT 


o o o cr c 

8.3856 


on c\o 

29. 9o 


MT-V 


HH 


8.2527 


31.347 




CHH 


8.2527 


31.357 




SVM 


8.2527 


31.360 




EIHH 


8.244 


31.358 




CRGC 




31.357 




FY 


8.2527 


31.364 




ATMS 




31.364 


MT-I/III 


HH 


8.5357 


30.310 




CHH 


8.5357 


30.31 




FY 


8.5357 


30.312 



these potentials are given as a sum of Gaussians but have a rather strong repulsion 
at short interparticle distances. This induces important two-body correlations in the 
wave function and consequently the convergence is slower than before. A grand angular 
momentum quantum number of 40 is sufficient for an accuracy of 1 keV in the binding 
energy, (iii) The MT-V potential is given as a superposition of Yukawians having a 
strong repulsive core with a 1/r divergence. This model, as well as the MT-I/III, 
represents the most challenging problem for the HH expansion, due to the difficulty of 
constructing accurate two-body correlations at short interparticle distances, where the 
cancellation between the kinetic and potential energy terms is critical. The introduction 
of correlation factors in the HH expansion allows for very accurate results in both the 
A = 3 and 4 systems. The HH without correlation factors requires a large number of 
basis elements to get convergence. For A = 3, the maximum value for the grand angular 
momentum quantum number (G max ) considered has been 200. For A = 4, it is possible 
to extrapolate the HH result to consider infinite number of basis elements, using the 
procedure of Ref. [65]. For the MT-V potential, the extrapolated HH result is 31.358 
MeV, very close to other accurate estimates. For the (s-wave) MT-I/III we observe that 
our estimate is already close to the very precise calculation of Ref. [99] . The "missing" 
binding energy in this case is estimated to be 21 keV, bringing our estimated binding 
energy to be 30.331 MeV. 
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6.2. Bound State with Realistic Interactions 

In this section, the A = 3 and 4 bound states are studied using the most recent models 
for the nuclear interaction, consisting of two- and three-nucleon potentials. We first 
summarize the main characteristics of these interaction models. 

Among the realistic models for the two-nucleon interaction, we have considered the 
Argonne AV8' [45], v 14 (AV14) [104] and AV18 [10] models, the Nijmegen 2 (NJ2) [9], 
and the CDBonn [11] models. We have also considered one of the models developed by 
Doleschall and collaborators [105], labelled as ISuj model [106]. Note that the Argonne 
and the NJ2 models are local and expressed in coordinate-space, the ISuj model is 
also expressed in coordinate-space but is non-local, and the CDBonn model is non- 
local and expressed in momentum-space. The AV8' potential is a reduction of the 
AV18 potential including only central, tensor and spin-orbit components, with modified 
parameters in order to reproduce the deuteron binding energy. We have included it, since 
it has been used before to produce a benchmark for the a-particle binding energy [66]. 
Among the potential models derived within an effective field theory approach, we have 
considered those ones calculated up to next-to-next-to-next-to-leading order by the 
Idaho [28] (N3LO-Idaho), and by the Jiilich [26, 27] (N3LO-Jiilich) groups. Note that the 
N3LO-Idaho model is obtained within a dimensional regularization scheme with cutoff 
parameter A=500 MeV. The N3LO- Jiilich model, instead, uses the spectral function 
regularization scheme with cutoff parameter A, and an exponential regulator function 
in order to remove the divergences in the Lippmann-Schwinger equation with cutoff 
parameter A. The cutoff combination [A, A] = [550, 600] MeV has been chosen in the 
present study. Finally, we have considered also two effective low-momentum interactions 
which have been derived from the AV18 model using a renormalization group approach. 
The smooth cutoff A; otu _fc of the exponential regulator function has been chosen equal 
to 2 or 3 fm _1 (Vi ow -k,2.o and Vi ow -k,3.o) [34]. Note that these two Vi ow -k models are 
different than the ones adopted in some of the earlier works on V[ ow -k [30, 31, 32, 33, 35], 
where a sharp cutoff formulation was used. 

Two remarks are here in order: (i) in the application of the HH expansion using 
the NJ2, ISuj, CDBonn, N3LO-Idaho, N3LO- Jiilich, and Vi ow _ k potential models, a 
truncation of the model space has been performed at j m ax=G, jmax being the maximum 
two-body total angular momentum considered. The chosen value allows for an accuracy 
of 1 keV for the triton binding energy, (ii) The point-like Coulomb interaction has been 
added when the 3 He and 4 He nuclei are considered. However, in the case of the AV18 
and NJ2 calculations, the full electromagnetic interaction has been added [9, 44]. 

As already discussed in Sec. 1, realistic models for the nuclear Hamiltonian include 
two-nucleon and three-nucleon interaction terms. Among the available two- and three- 
nucleon interaction models, we have used in the present work the AV18 together with 
the Urbana IX [14] (AV18/UIX) model and the CDBonn together with the Tucson- 
Melbourne [16] model (CDBonn/TM). In the TM model, the regularization cutoff 
parameter Atm has been fixed to 4.795 m n , in order to reproduce the triton binding 
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Table 5. The triton and the 4 He binding energies B (MeV), calculated with the AV18, 
CDBonn, N3LO-Idaho, two-nucleon interaction models, and with the AV18/UIX, 
CDBonn/TM, and N3LO-Idaho/N2LO two- and three-nucleon interactions. Note that 
in the CDBonn/TM case, the parameter Atm has been chosen to be 4.795 m n for 3 H 
and 4.784 m T for 4 He. The HH results are compared with the ones obtained with 
different approaches. See text for explanations. Note that a recently updated FE/FY 
result [112] for the 4 He binding energy calculated with the CDBonn potential is given 
in parentheses. 



Potential 


Method 


B( 3 H) 


5( 4 He) 


AV18 


PHH 


7.624 


_ 




HH 


7.624 


24.22 




FE/FY [108, 107] 


7.621 


24.23 




FE/FY [110] 


7.621 


24.24 




FE/FY [39] 


7.616 


24.22 


CDBonn 


HH 


7.998 


26.13 




FE/FY [108, 107, 112] 


8.005 


26.23 (26.16) 




FE/FY [111] 


7.998 


26.11 




NCSM [56] 


7.99(1) 




N3LO-Idaho 


HH 


7.854 


25.38 




FE/FY [112] 


7.854 


25.37 




FE/FY [111] 


7.854 


25.38 




NCSM [109] 


7.852(5) 


25.39(1) 


AV18/UIX 


PHH 


8.479 






HH 


8.479 


28.47 




FE/FY [108, 107] 


8.476 


28.53 




FE/FY [39] 


8.473 




CDBonn/TM 


HH 


8.474 


29.00 




FE/FY [108, 107] 


8.482 


29.09 


N3LO-Idaho/N2LO 


HH 


8.474 


28.37 




NCSM [56] 


8.473(5) 


28.34(2) 



energy when used in conjunction with the CDBonn potential. The parameters a, 6, c 
and d of the potentials are the ones of Table XI of Ref. [107]. However, in the 4 He 
calculation, for the sake of comparison, the CDBonn/TM model has been used with 
Atm = 4.784 m n , the same as in Ref. [108]. 

Within an effective field theory approach, three-nucleon interactions arise when 
next-to-next-to leading order (N2LO) contributions are considered. Therefore, the 
A = 3 and 4 nuclei have been studied also using the N3LO-Idaho two-nucleon together 
with the local N2LO three-nucleon interaction as developed by Navratil in Ref. [109] 
(N3LO-Idaho/N2LO). The two free parameters in the N2LO three-nucleon interaction 
model have been chosen from the combination that reproduces the A = 3, 4 binding 
energies [109]. 

The A = 3 and 4 binding energies and ground-state properties are listed in Tables 5, 
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Table 6. The triton binding energies B (MeV), the proton r p and neutron r„ 
radii (fm), the expectation values of the kinetic energy operator (T) (MeV), the 
mixed-symmetry S', P, D and T = 3/2 probabilities (all in %) calculated with the 
AV8', AV14, AV18, NJ2, ISuj, CDBonn, N3LO-Idaho, N3LO-Jiilich, V low - ki2 .o and 
Vi ow -k,3.o two-nucleon interaction models, and with the AV18/UIX, CDBonn/TM 
(with Atm = 4.795 and N3LO-Idaho/N2LO two- and three-nucleon interactions. 



Potential 


B 


(T) 


r P 


r n 


Ps> 


Pp 


Pd 


Pt=3/2 


AV8' 


7.767 


47.605 


1.642 


1.807 


1.273 


0.067 


8.579 


0.0 


AV14 


T O O A 

7.684 


45.678 


1.667 


1.828 


1 1 or* 

1.126 


0.076 


8.967 


0.0 


AV18 


7.624 


46.727 


1.653 


1.824 


1.293 


0.066 


8.510 


0.0025 


NJ2 


7.651 


47.520 


1.647 


1.816 


1.292 


0.064 


8.329 


0.0032 


ISuj 


8.475 


32.950 


1.568 


1.715 


1.433 


0.026 


4.776 


0.0028 


CDBonn 


7.998 


37.630 


1.618 


1.771 


1.310 


0.047 


7.018 


0.0049 


N3LO-Idaho 


7.854 


34.555 


1.655 


1.808 


1.365 


0.037 


6.312 


0.0009 


N3LO-Jiilich 


7.292 


43.467 


1.726 


1.892 


1.690 


0.018 


4.316 


0.0017 




8.597 


29.509 


1.574 


1.713 


1.330 


0.018 


4.062 


0.0016 




8.085 


32.790 


1.614 


1.765 


1.317 


0.042 


6.749 


0.0020 


AV18/UIX 


8.479 


51.275 


1.582 


1.732 


1.054 


0.135 


9.301 


0.0025 


CDBonn/TM 


8.474 


39.364 


1.580 


1.722 


1.202 


0.101 


6.971 


0.0049 


N3LO-Idaho/N2LO 


8.474 


36.482 


1.611 


1.752 


1.242 


0.121 


6.815 


0.0009 


Exp. 


8.482 




1.60 













6, 7, and 8. In Table 5, using some representative models among the ones quoted 
above, i.e. the AV18, CDBonn, N3LO-Idaho, AV18/UIX, CDBonn/TM and N3LO- 
Idaho/N2LO, the PHH and HH results for the triton and 4 He binding energies are 
compared with those ones obtained with other approaches. In particular, we have 
referred to the FE (for A = 3) and FY (for A = 4) approach of Refs. [108, 110, 111, 39], 
and the NCSM approach of Ref. [54, 56, 109]. 

From inspection of Table 5, we can conclude that for any of the considered realistic 
potential models, the agreement between the HH or PHH results and those from the 
other techniques here considered is excellent. Such an agreement is present for both 
local and non-local potentials. To be noticed that the AV18 result of Ref. [39] does not 
include the T = 3/2 contribution. Moreover, the CDBonn and CDBonn/TM results of 
Ref. [108] include the n — p mass difference. Once this contribution is subtracted, the 
agreement with the other methods is quite good. In the same row of the table a recently 
updated result for 4 He is given in parentheses [112]. 

In Tables 6, 7, and 8 we present our results for the A = 3, 4 bound states properties, 
among which the expectation value of the kinetic energy, the proton and neutron radii 
and different occupation probabilities. In the case of inclusion of NN forces only, 
we notice that all the potentials, except for the ISuj model, do not reproduce the 
experimental trinucleon binding energy. The ISuj model has been constructed modifying 
the off-energy shell part of AV18 in order to reproduce the NN scattering data and also 
some A = 3 observables, as the trinucleon binding energy. This model predicts a quite 
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Table 7. Same as Table 6, but for 3 He. 



Potential 


B 


(T) 


r p 


r n 


Ps> 


Pp 


Pd 


Pt=3/2 


AV8' 


7.108 


46.696 


1.846 


1.666 


1.470 


0.066 


8.542 


0.0 


AV14 


7.033 


44.813 


1.867 


1.691 


1.315 


0.076 


8.967 


0.0 


AV18 


6.925 


45.685 


1.872 


1.678 


1.530 


0.065 


8.467 


0.0080 


NJ2 


6.994 


46.607 


1.857 


1.668 


1.491 


0.063 


8.282 


0.0083 


ISuj 


7.712 


32.194 


1.757 


1.590 


1.643 


0.026 


4.795 


0.0059 


CDBonn 


7.263 


36.767 


1.819 


1.637 


1.542 


0.046 


7.000 


0.0109 


N3LO-Idaho 


7.128 


33.789 


1.855 


1.675 


1.607 


0.037 


6.313 


0.0062 


N3LO-Jiilich 


6.590 


42.471 


1.948 


1.751 


1.985 


0.018 


4.327 


0.0095 




7.847 


28.861 


1.755 


1.591 


1.539 


0.018 


4.086 


0.0063 




7.358 


32.087 


1.810 


1.632 


1.537 


0.042 


6.736 


0.0073 


AV18/UIX 


7.750 


50.211 


1.771 


1.602 


1.242 


0.132 


9.248 


0.0075 


CDBonn/TM 


7.720 


38.495 


1.767 


1.597 


1.409 


0.099 


6.966 


0.0106 


N3LO-Idaho/N2LO 


7.733 


35.745 


1.794 


1.628 


1.450 


0.119 


6.818 


0.0057 


Exp. 


7.718 




1.77 













low D-state probability as the models derived from %PT. The binding energy and the 
other quantities predicted by the Vi ow -k models depend sizably on Ai ow -k- Here, we 
have chosen two representative values of such a parameter. 

In the cases in which a 3N interaction is included, the experimental binding energy 
is well reproduced, since one of the parameter of such a term is chosen accordingly. In 
all the calculations the n — p mass difference contribution has not been included. This 
contribution can be estimated perturbatively by taking the mean value of the operator 

^ = Ej (4 " i) v? T « (i) ' (118) 

This effect has been found to be quite tiny. For example, using the AV18 interaction 
model, the change in binding energy of 3 H, 3 He, and 4 He has been found to be about +6 
keV, —6 keV, and —0.15 keV, respectively [108, 65]. By inspecting the tables, it can be 
observed the presence of quite large differences in the occupation probabilities and the 
expectation value of the kinetic energy for the different models. For example, the mean 
value of the kinetic energy for CDBonn, ISuj, N3LO-Idaho and Vi ow -k-type potentials 
is noticeably smaller than that found with the AV18 or NJ2 potentials. This is due 
to the fact that the repulsion at short interparticle distances is softer for the former 
potentials than for the latter ones. Also the percentages of the P- and .D-waves are 
significantly smaller for the former potentials. This fact has an important consequence 
in the application of the HH method, since for potentials with a softer repulsion at short 
interparticle distances the convergence is usually faster. For example, for A = 4, the 
convergence has been reached including HH states up to G max = 80, 52, 32, 22 for AV18, 
CDBonn, N3LO-Idaho, and Vi ow -k potential models, respectively. 

The protonic radii are well reproduced by the models including a 3N interaction. To 
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be noticed that the ISuj and Vi ow _k,2.o models, which reproduce quite well the A = 3, 4 
binding energies with only two-body forces, tend to predict too small protonic radius. 
This could be related to the corresponding small D-wave percentage. 

A few comments on the estimated percentages of the T = 3/2 and T > 
components in the A = 3 and 4 systems are in order. For the percentage of the T = 1 
component in the 4 He ground state wave function, 50% of it is due to the effect of 
the Coulomb interaction between the protons, the remaining 50% is due to the CSB 
terms in the nuclear interaction. This can be seen from Table 8, since the models 
where the CSB is absent (AV8' and AV14) predict a percentage approximately half of 
those calculated with AV18, CDBonn, N3LO-Idaho, etc., which were fitted to both pp 
and np data including CSB by construction. The difference in Pt=i for these latter 
models is however small, due to the important role of the Coulomb interaction. On 
the contrary, the T = 3/2 component in triton and the T = 2 component in 4 He are 
largely dominated by CSB of nuclear origin (different pion masses, etc.). The values 
reported in Table 6 and 8 show that, depending on the interaction, rather different 
values for the triton Pt=3/2 and the 4 He Pt=2 are obtained (note that the standard 
models of 3N have little effect on the isospin admixtures [65]). The origin of the rather 
large differences found for the triton Pt=3/2 and the 4 He Pt=2 (a factor 5 between 
CDBonn and N3LO-Idaho) must be related to quite different off-shell behavior of the 
CSB terms of the interactions. The T = 3/2 component in 3 He is also affected by the 
Coulomb interaction, which strongly reduces the difference between the Pt=3/2 results. 
The knowledge of the T = 1 and 2 percentages is important for parity-violating electron 
scattering experiments on 4 He, aimed at studying admixture of strange quark ss pairs 
in nucleons and nuclei [113, 114, 115]. A preliminary study of this important aspect has 
been published in Ref. [116]. It could play an important role also in the study of the 
reaction d + d — > a + 7r°. This reaction is possible only if isospin symmetry is violated, 
namely it probes directly the CSB terms in the nuclear Hamiltonian [117, 118]. 

Finally, we would like to comment on the capability of these potential models in 
reproducing simultaneously the A = 3, 4 binding energies. It is well known that there 
is an almost linear relation between the 3 H and 4 He binding energies. By inspection of 
Tables 6, 7, and 8 we observe that the models which reproduce the 3 H binding energy 
(ISuj, AV18/UIX, CDBonn/TM, and N3LO-Idaho/N2LO) slightly overpredict the 4 He 
binding energy. This overprediction is of about 0.6 MeV for ISuj and CDBonn/TM and 
reduces to 0.16 MeV for AV18/UIX or to 0.06 MeV for N3LO-Idaho/N2LO. Another 
important aspect is the correct description of the mass difference between the 3 H 
and 3 He, directly related to the CSB components in the Hamiltonian. As already 
mentioned, the present calculations do not include the n — p mass difference. Taking 
it into account together with other small contributions, as explained in Ref. [108], the 
AV18/UIX theoretical estimate for the mass difference D is around 750 keV, slightly 
smaller than the experimental one D cxp = 764 keV. The 3 H- 3 He mass difference values 
for the CDBonn/TM and N3LO-Idaho/N2LO models are D = 754 keV and 741 keV, 
respectively, as can be seen from Tables 6 and 7. Taking into accout the n — p mass 
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Table 8. The 4 He binding energies B (MeV), the proton r p and neutron r„ radii 
(fm), the expectation values of the kinetic energy operator (T) (MeV), the P, D, 
T = 1, and T = 2 probabilities (all in %) calculated with the AV8', AV14, AV18, 
NJ2, ISuj, CDBonn, N3LO-Idaho, N3LO-Julich, Vi ow -k,2.o and Vi ow -k,3.o two-nucleon 
interaction models, and with the AV18/UIX, CDBonn/TM (with A T m = 4.784 m,), 
and N3LO-Idaho/N2LO two- and three- nucleon interactions. 



Potential 


B 


(T) 


r p 


r„ 


Pp 


Pd 


Pt=i 


Pt=2 


AV8' 


24.90 


98.80 


1.507 


1.502 


0.408 


14.50 


0.0013 


0.0001 


AV14 


24.23 


95.76 


1.526 


1.521 


0.393 


14.41 


0.0013 


0.0001 


AV18 


24.21 


97.84 


1.514 


1.509 


0.347 


13.74 


0.0028 


0.0052 


NJ2 


24.42 


100.27 


1.506 


1.501 


0.334 


13.37 


0.0016 


0.0074 


ISuj 


28.97 


68.66 


1.383 


1.378 


0.097 


6.10 


0.0020 


0.0072 


CDBonn 


26.13 


77.58 


1.458 


1.453 


0.223 


10.74 


0.0029 


0.0108 


N3LO-Idaho 


25.38 


69.24 


1.518 


1.513 


0.172 


9.29 


0.0035 


0.0024 


N3LO-Jiilich 


22.20 


86.85 


1.636 


1.631 


0.081 


6.22 


0.0065 


0.0063 




29.66 


60.77 


1.389 


1.384 


0.064 


4.69 


0.0020 


0.0023 




26.57 


65.68 


1.455 


1.450 


0.201 


10.04 


0.0025 


0.0038 


AV18/UIX 


28.46 


113.30 


1.430 


1.425 


0.732 


16.03 


0.0025 


0.0050 


CDBonn/TM 


29.00 


84.56 


1.396 


1.391 


0.454 


9.94 


0.0021 


0.0105 


N3LO-Idaho/N2LO 


28.36 


74.93 


1.476 


1.471 


0.608 


10.79 


0.0028 


0.0020 


Exp. 


28.30 




1.47 













difference contribution 12 keV), it can be seen that the CDBonn/TM well reproduces 
.Dexp, whereas the N3LO-Idaho/N2LO model gives a slightly smaller value. For the ISuj 
model, the value D = 763 keV coincides with _D exp . To be noticed that additional 
contributions could come from CSB terms in the 3N interaction. However, the 3N 
models considered in the present analysis do not contain terms which violate isospin 
symmetry. 

6.3. A = 3 and 4 Scattering States 

The A = 3 and 4 zero-energy scattering states have been studied with the following 
interaction models: the two-nucleon AV14, AV18 and N3LO-Idaho, and the two- and 
three-nucleon AV18/UIX and N3LO-Idaho/N2LO models. In a consistent calculation of 
the A = 3,4 scattering lengths, the electromagnetic interaction cannot be disregarded. 
In particular, the magnetic moment interaction gives contribution in processes in which 
the incident nucleon is either a proton or a neutron. In the first case, besides the 
Coulomb and the magnetic moment interactions, other contributions of long range as 
vacuum polarization and two-photon-exchange are present. The treatment of these 
terms remains outside the aims of the present work. Accordingly, in the description of 
pd and p 3 He, we consider only the point Coulomb interaction and, in some cases, the 
magnetic moment interaction too. 

The results for the nd and pd doublet and quartet scattering lengths are given in 
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Tables 9 and 10 and are compared with those obtained by other techniques [119, 120] 
and the available experimental data [121, 122]. The pd calculation has been performed 
with and without the inclusion of the T = 3/2 components of the wave function 
(see Eq. (102)). Furthermore, we performed two different calculations using the AV18 
potential since, by construction, it includes the long range electromagnetic terms. In 
the nd case, the results labelled as AV18* do not include these terms. In the pd case 
the AV18* results include only the point Coulomb interaction, whereas those labelled as 
AV18 include the complete electromagnetic potential with the exception of the two- 
photon exchange and the vacuum polarization terms. The calculations have been 
performed with the PHH expansion in the A = 3 case with local potential models, 
and with the HH expansion in all other cases. The dimension of the basis and the 
number of the integration points have been increased in order to reach an accuracy of 
0.002 fm in the calculation of the scattering lengths. 

From inspection of the tables we can conclude that: (i) both the nd and pd 
quartet scattering lengths are almost model-independent, (ii) The nd doublet scattering 
length 2 a n d is very sensitive to the choice of the two-nucleon potential model, when 
no 3N interaction is included. Once the 3N interaction is included, and therefore 
the triton binding energy is well reproduced, 2 a n d becomes little model-dependent. 
This is a well-known feature, related to the fact that 2 a n d and the triton binding 
energy are linearly correlated (the so-called Phillips line [123]). (iii) The pd doublet 
scattering length 2 a p d is positive and quite model-dependent, if only the two-nucleon 
interaction is included. Once the three-nucleon interaction is added, 2 a p d becomes 
small. Some model-dependence remains, but the problem of extrapolating to zero 
energy the experimental results makes very difficult any meaningful comparison between 
theory and experiment, (iv) The nd quartet scattering length A a n d is in very good 
agreement with the experimental data. Conversely, some disagreement is observed in 
the doublet scattering length 2 a nd for most of the NN plus 3N interactions considered. 
The results presented have been calculated without including the neutron-proton mass 
difference. This contribution will further decrease the value of 2 a n d- To be noticed 
that a recent measurement [122] of the coherent nd scattering length predicts 2 a n d = 
[0.645 ± 0.003(expt) ± 0.007(theory)] fm. This very precise measurement is not well 
described by any of the potential models analyzed in this section, with the exception of 
the N3LO-Idaho/N2LO. 

The results for the n 3 H and p 3 He singlet and triplet scattering lengths are given in 
Table 11 and are compared with other theoretical calculations, performed using the FY 
equations method [39, 111], and the available experimental data [124, 125, 126, 127]. 
For n 3 H, the scattering lengths can be obtained from the experimental values of the 
total cross section and the coherent scattering length a c , 



The n 3 H cross section has been accurately measured over a wide energy range and 
the extrapolation to zero energy does not present any problem. The value obtained is 




a, 




a . 



(119) 
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Table 9. The nd doublet and quartet scattering lengths 2 a n d 1 4 a nd , in fm, are 
calculated with the AV14, AV18, AV18*, and N3LO-Idaho two-nuclcon interaction 
models, and with the AV14/TM (with A TM = 5.13 m n ), AV18/UIX, AV18*/UIX, 
and N3LO-Idaho/N2LO two- and three-nucleon interactions, and are compared with 
the results obtained with the FE technique [119, 120] and the available experimental 
data [121, 122]. 



Potential 


Method 


2 


4 

a-nd 


AV14 


PHH 


1.189 


6.379 


AV14 


FE [119] 


1.204 


6.380 


AV18 


PHH 


1.258 


6.345 


AV18 


FE [120] 


1.248 


6.346 


AV18* 


PHH 


1.275 


6.325 


AV18* 


FE [120] 


1.263 


6.326 


N3LO-Idaho 


HH 


1.100 


6.342 


AV14/TM 


PHH 


0.586 


6.371 


AV18/UIX 


PHH 


0.590 


6.343 


AV18/UIX 


FE [120] 


0.578 


6.347 


AV18*/UIX 


PHH 


0.610 


6.323 


AV18*/UIX 


FE [120] 


0.597 


6.326 


N3LO-Idaho/N2LO 


HH 


0.675 


6.342 


Exp. [121] 
Exp. [122] 


0.65±0.04 6.35±0.02 
0.645±0.003±0.007 



Table 10. Same as Table 9 but for pd. Also listed are the results obtained without 
the inclusion of the T — 3/2 component in the wave functions ( 2 a n d(T = 1/2) and 
i a nd {T= 1/2)). 



Potential 


Method 


2 a P d 


4 

CLpd 


2 a pd (T = 1/2) 


4 Opd(T = 1/2) 


AV14 


PHH 


0.937 


13.773 


0.941 


13.773 


AV14 


FE [119] 






0.965 


13.764 


AV18 


PHH 


1.134 


13.662 


1.150 


13.662 


AV18* 


PHH 


1.185 


13.588 


1.198 


13.589 


N3LO-Idaho 


HH 


0.876 


13.646 


0.866 


13.646 


AV18/UIX 


PHH 


-0.089 


13.662 


-0.074 


13.663 


AV18*/UIX 


PHH 


-0.035 


13.588 


-0.019 


13.590 


N3LO-Idaho/N2LO 


HH 


0.072 


13.647 


0.082 


13.647 



o~t = [1.70±0.03] b [128]. The coherent scattering length has been measured by neutron- 
interferometry techniques. The most recent values reported in the literature are a c = 
[3.59 ± 0.02] fm [124] and a c = [3.607 ± 0.017] fm, the latter value being obtained from 
p 3 He data using an approximate Coulomb-corrected R-matrix theory [125]. However, in 
the 1 a- 3 a plane, the ellipse defined by the total cross section and corresponding to the 
experimental value of o T = 1.70 b and the straight line corresponding to the coherent 
scattering length a c 3.6 fm are almost tangent. Therefore, a small change in the 
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a c value produces a large variation of 1 a and 3 a. This is also the reason for the large 
uncertainty in the values reported in Table 11. 

Note that there is a fairly good agreement between the results of the various 
theoretical calculations. The calculated scattering lengths show a scaling behaviour with 
respect to the 3 H binding energy similar to the Phillips line, as discussed in Ref. [72]. The 
total cross sections obtained with the AV18/UIX and N3LO-Idaho/N2LO models are 
1.73 b and 1.68 b, in agreement with the experimental value. The theoretical coherent 
scattering lengths are 3.73 fm and 3.65 fm, respectively, at variance with respect to the 
experimental value. The origin of this discrepancy is still unclear. 

For p 3 He, scattering length measurements were reported in Refs. [126, 127]. 
However, these have rather large errors. In fact, the p 3 He data have been extrapolated 
to zero energy from measurements taken above 1 MeV, and therefore suffer of large 
systematic uncertainties. 

Table 11. The n 3 H, p 3 He singlet and triplet scattering lengths 1 a n 3 H , 3 a„3 H , 1 a p 3 He , 
3 a p 3 Ho , in fm, are calculated with the AV18, and N3LO-Idaho two-nucleon interaction 
models, and with the AV18/UIX and N3LO-Idaho/N2LO two- and three- nucleon 
interactions. The available results from the FY equations method [39, 111] and the 
available experimental data [124, 125, 126, 127] are also reported. 



Potential 


Method 


a n 3 H 


3 a„3 H 


1 a p 3 Hc 


3 


AV18 


HH 


4.29 


3.73 


12.9 


10.0 


AV18 


FY [39] 


4.27 


3.71 






AV18 


FY [111] 


4.28 


3.71 






N3LO-Idaho 


HH 


4.20 


3.67 


11.5 


9.2 


N3LO-Idaho 


FY [111] 


4.23 


3.67 






AV18/UIX 


HH 


4.10 


3.61 


11.5 


9.1 


AV18/UIX 


FY [39] 


4.04 


3.60 






N3LO-Idaho/N2LO 


HH 


3.99 


3.54 


11.0 


8.6 


Exp. [124] 




4.98 ±0.29 


3.13 ±0.11 






Exp. [125] 




4.45 ±0.10 


3.32 ±0.02 






Exp. [126] 








10.8 ±2.6 


8.1 ±0.5 


Exp. [127] 










10.2 ± 1.5 



7. Summary and Conclusions 

In this article we have presented a theoretical approach suitable to describe with high 
precision bound states and zero-energy scattering states for the A = 3, 4 nucleon 
systems. A particular attention has been given to the presentation of the various models 
available for the two- and three-body nuclear interaction. In fact, in recent times, 
different NN and 3N interactions have been proposed and applied by many authors to 
the study of few-nucleon systems. The capability of these interactions to describe all the 
complexity of the nuclear dynamics is at present an intense subject of research. Bound 
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states and low-energy scattering states are the first structures to be examined using the 
interaction models under consideration. To this aim, any reliable technique employed 
should allow for meaningful comparison with the available experimental data. In other 
words, the solution of the Schrodinger equation should be accurate enough in order to 
eliminate uncertainties in the theoretical predictions. 

The method discussed in this paper is based on the expansion of the wave function 
over the HH functions and the solution of the Schrodinger equation is obtained by using 
a variational principle. Nowadays computing facilities make it possible to consider 
very large sets of HH functions, so that the problem of convergence can be carefully 
studied for specific observables as binding energies or scattering lengths. Several 
papers [61, 62, 65, 67, 68] have been devoted to applying and verifying this approach 
for different systems, and the obtained results have been successfully compared to those 
of other accurate techniques. 

There are two main motivations behind the present work. The first one is to 
highlight the details of the HH method and its capability and accuracy in studying 
A = 3 and 4 nuclear systems. The extension to larger A-values is an important problem 
and studies in this direction are currently in progress. The second motivation is that 
the results presented here for different types of interaction models should turn out to 
provide a useful comparison for other available or future techniques. The main merit 
of the HH method is that it can be implemented to treat the A = 3,4 systems with 
interactions given in configuration or in momentum space including the electromagnetic 
potential. In particular, the description of low-energy scattering states with charged 
particles does not arise particular difficulties and, in fact, the HH technique has been 
applied successfully for the description of the pd capture reaction [77, 80] and the p 3 He 
weak capture reaction [78, 79] in the keV region. 

From the results here presented we can conclude that there are some problems, 
not yet completely solved, in the theoretical description of the A = 3, 4 systems. 
The simultaneous reproduction of the 3 H, 3 He and 4 He binding energies is not totally 
satisfactory. Several combinations of modern NN plus 3N interaction models tend to 
overestimate the a-particle binding energy, once the 3N interaction strength has been 
fixed to reproduce the triton binding energy. The nucleons in the a-particle are in 
average very close to each other, so that a greater sensitivity to 3N interaction terms 
should be manifested in this system. The N2LO 3N interaction model in conjunction 
with the N3LO-Idaho NN interaction gives the better description of the three binding 
energies mentioned above. This model contains a number of operators very similar to 
those present in the TM and UIX model, but having different relative strengths. The 
overestimation given by the CDBonn/TM potential is appreciable, suggesting that a 
further analysis of the TM model could be useful. 

Another problem is the underestimation of the nd doublet scattering length by all 
the models considered except the N3LO-Idaho/N2LO. A high precision measurement of 
the coherent nd scattering length is available and, from this quantity, the nd doublet 
scattering length can be evaluated. This has been done using the theoretical estimates 
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for the nd quartet scattering length since this quantity results to be almost model 
independent. In any case, after considering the theoretical uncertainty introduced, the 
theoretical estimates given in Table 9 are lower with respect to the mean experimental 
value by 5 % to 10 %. This is not the case for the N3LO-Idaho/N2LO potential. 
It seems that this model is the only one among those studied here able to reproduce 
simultaneously the experimental values of the A = 3 and 4 binding energies and the 
A = 3 scattering lengths. 

Although many efforts have been done, and are still in progress, for a more accurate 
determination of the nuclear interaction, a long way seems to be in front of us. To this 
aim, the extension of the HH method to treat systems with A > 4 and scattering states 
at medium and high energies will be very useful. Intense efforts are pursued in this 
directions by the authors. 



Appendix A. The Transformation Coefficients for A = 3 

Let us discuss how to compute the transformation coefficients (TC) for A = 3. It is 
convenient to start from Eq. (70), which can be written as 

"Sana^nJ = / ^2 P^J^J.LM^)] 1 ^H^n^^M^) , (A.l) 

where the functions ^H{e 1 e 2n2 },LM( l ^2) are defined in Eq. (53). First of all, let us recall 
the following identity satisfied by the spherical harmonics 

c % m (c)= E «H f ^% i4 fy,(fl)y 4 (6)| , (A.2) 

where the vector variables are related by c = a + b and 



(21+1)! rA ,x 
ViMb ~ " (2£ a + l)!(24+l)! ' [ ' 

The Jacobi vectors constructed with permutation p are linearly related to the Jacobi 
vectors xi, x 2 for the permutation 1 (corresponding to the order 1, 2, 3 of the particles), 
namely 

M ~ E > * = 1 > 2 > (A-4) 



.7=1,2 



where otf) are numerical coefficients. It can be shown that 
(sin^VHcos^V 2 W{xf)Y i2 {xf) = 

L J L,M 

E C $2,*i* 2 ( Sin ^ 2 ' COS ^ 2 ) [ y 4(^l)^a(^2)] L)M , ( A - 5 ) 

where 
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x D tlMM D ia ^(-) Xl+ ^ + ^iJJ' 1 ^X 1 X 2 x: i ^ 

■(iS5)(i*i){iiiH 

and 1 = V21TT. 

The integrand in Eq. (A.l) depends on the angular variables Xi and x 2 through 
the spherical harmonics and the argument cos 2ip y f > of the Jacobi polynomial, which is 
given by 

\a$)\ f m<p 2 ) 2 + (a$) 2 (co a <p 2 ) 2 

1 - 1 , (A.7) 



cos 2ip^ 



(v) (p) 

+ 2a 21 ct 22 fJ> sin ip 2 cos ip 2 



where \i = x-y • x 2 . It is possible to perform analytically part of the angular integration 
(keeping fx fixed) with the result 

J d Xl dx 2 [^(xi)^^)]^ [^ 1 (xi)y /a (x 2 )] £)M = 



where 



Af 2 ^ 2 ' L = (-) L+ * 2+ ^f/ 2 f/ 2 (2A + 1) 

' t to. A 






Finally, using Eqs. (A. 6) and (A.8), we arrive at the expression for the TC 



(A.8) 



xP n+V2^ + i/2 (cos2 ^ )p ^ +1 / 2A+1/2(cos2 ^)^ 



qf4 iAiA2 (sinv. 2 ,cosv. 2 )Ar 2A1A2 "PA(/i) , (A.9) 

Ai,A2,A 

where A/^ 2 '^ 2 is given in Eq. (34) and v 2 = G + 2. This integral can be calculated easily 
using Gauss quadrature, since it is a polynomial in fx, cosy? 2 and sin</? 2 . 



Appendix B. The Transformation Coefficients for A = 4 

Let us now consider the calculation of the transformation coefficients (TC) for A = 4 
(N = 3) for a generic permutation p of the particles. In the following, the hypervariables 
constructed starting from a generic choice of the Jacobi vectors x^, x 2 p \ x^, 
corresponding to the permutation p of the particles and either set A or B (see Eq. (6)), 
will be denoted by We will also consider the hyperangular variables f2 3 constructed 
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in terms of a reference set of Jacobi vectors x 1; x 2 , x 3 corresponding for example to those 
defined in set A and to the order 1, 2, 3, 4 of the particles. Let us define 

(B.l) 



y P 



cos 2<4 p) 



cos 20^ 



where the hyperangles <j) 2 



(p) 



are defined in Eq. (8) in terms of the moduli of the 



Jacobi vectors, 
function is 



In terms of these variables, the expression of a generic A = 4 HH 



ri {£ 1 ,£ 2 ,e 3 ,L 2 ,n 2 ,n 3 },LM\ iL 3 ) 



(Y ei (x?)Y £2 (x?)) L Y e3 (x 



(ph 

3 i 



LM 



X M h/2 ' 
~ Jv n 2 ,n 3 



(1-^)^(1+^)" 



X (1 



2 (1 + Z V 



x R 



712 



nz 



(z p ) , (B.2) 



where the normalization factor is 

\ ^i+^2+n2+^- 



Jx n 2 ,n 3 I ~ 



X 



n 



i=2 



r(y. - n 3 - tj - l/2)r(n, + i 3 + 3/2) J ' (R3) 

with Vj = Gj + (3j - 5)/2 and Gj defined in Eq. (30). The superscript G in Eq. (B.2) 
has been inserted to remember that the possible choices of the quantum numbers 
{h, ^2, h, L 2 ,n 2 ,n 3 } are restricted to the case t\ + £ 2 + £3 + 2(n 2 + n 3 ) = G. 

Let us start with the state having quantum numbers n 2 = n 3 = (G = G m = 
£1 + £2 + £3)- Since the values of G m to be considered are small, the TC for this state 
can be easily calculated by means of any of the methods proposed in the literature. For 
example, one can employ the technique of Ref. [129]. This method is based on the fact 
that Eq. (69) has to be verified for any spatial configuration of the Jacobi vectors set 
{xj}. For given values of G and L, the sum in the right-hand side of Eq. (69) is over 
N GL terms. If the equation is required to be satisfied for N GL values of the set {x;}, 
one gets a system of N GL linear equation for the required TC. The HH functions we 
are interested in have G = G rn = £\ + 1 2 + £ 3 < 6, therefore the number N G l of linear 
equations is small and a standard numerical technique can be employed. 

The Jacobi vectors x^, x^, x^ are linearly related to the Jacobi vectors xi, x 2 , x 3 , 
hence the following relations hold 



( x 3 



(p) 



f = E rj 



(p) 



Xj • X 



3 1 



(4 P) ) 2 = E A 



ij ' 



(B.4) 



*J=1 



where the numerical coefficients and can be easily computed by expressing the 
vectors xff"* and x^ in terms of Xi, x 2 and x 3 . 



» 



Let us now assume to know the TC for the HH functions ^Tif ( 
in terms of the HH functions constructed with f2 3 , namely 



{£i,£ 2 ,£ 3 ,L 2 ,n 2 ,n 3 },LM 



(4) W 



G 

{e 1 ,e 2 ,£ 3 ,L 2 ,n 2 ,n 3 },LM 



E 



a 



(p),G,L 

£\ ,£ 2 ,£ 3 ,L 2 ,n 2 ,n 3 \£\ ,£' 2 ,£' 3 ,L' 2 ,n' 2 ,n' 3 
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X (4) ^' 1 ,^,4,L' 2 ,^^},LM( n 3) • (B.5) 

It is important to note that the sum over {£[, £' 2 , £' 3 , L 2 , n' 2) n' 3 } is restricted by the 
condition £[ + £' 2 + £' 3 + 2(n' 2 + n 3 ) = G. This can be understood from the fact that 
p G x ^7i G is a harmonic polynomial of degree G. The transformation given in Eq. (B.4) 
clearly cannot change the degree of such a polynomial and therefore ^T-C G (Q^) can be 
expressed only in terms of HH functions with the same G. Also the transformation (B.4) 
is equivalent to a rotation, and therefore L, M cannot change and the TC cannot depend 
on M. 

Let us consider first the function (4) ^L^,i 2 ,n2,n 3 +i},LM( fi 3 P) )- Using the 
expression of the Jacobi polynomials P n . 3+ i(z p ) in terms of P n3 (z p ) and P n . 3 ^i(z p ), one 
obtains 

(p),G+2,L 
u t 1 ,£ 2 /3 M ,n 2 ,n 3 +l;e , 1 ,t' 2 /' 3 ,L 2 ,n' 2 ,n' 3 

3 

2b n3 ^ a i 1 } 2 IhM^n^i'l^t^L'^n'l XI r *>J ^'3 ' ( B ' 6 ) 

^1 '^2 '^3 1 ^'2 '^3 ^'J ^" 

where 

t _ ■^n 2 ,n 3 +l 

(2n3 + a3 + /33 + l)(2n 3 + a3 + /3 3 + 2) 
3 ~ Milfnz 3 2(n 3 + l)(n 3 + a 3 + ft + 1) ' 1 J 

and a 3 = £i + £ 2 + 2n 2 + 2, j3 3 = £ 3 + 1/2. The terms ijj, i, j = 1, 3 are defined by 

t x, ; • x, 



where dQ 3 is given by 



P 2 

X (4) ^§iA/ 3 ,L2,n 2 ,n 3 },iM( fi 3) , (B.8) 



dtt 3 = — X —=dx 1 dx 2 dxj,dydz(\ - j/)3(l + - ^) 2 (1 + z)s , (B.9) 

128V 2 

and y = cos20 2 , z = cos20 3 . The integrals Iij involve only functions constructed within 
the same set of Jacobi vectors Xj, % — 1, 3. Moreover, the factors (xj • Xj)/p 2 , have 
the following expressions in terms of the hyperangular variables f2 3 = (x±, x 2 , £3, y, z) 
(remember that x 3 = pcosip 3 , etc), 

v 11 x? (l-y)(l-*) r . 

i,J = M, ^ = ^ , (B.10) 

99 ^ (l + y)(l-*) 
i,J = 2,2, = , (B.ll) 

U = 3,3, ij = 1 ^> (B-12) 
p z 2 



. . xi • x 2 4vr (1 - z)\Jl — y 2 , . 

= 1.2 , — — = --^= [yi(xi)yi(a; 2 )] 0)0 , (B.13) 



•■ 10 X 1' X 3 47T VI - ^V 1 -V r W ~ ^ v( * \\ , R1 ,N 

= 1,3 , — ^— = --^= ^= [Yi(xi)yi(x 3 )] 00 , (B.14) 

• ■ oq X2 ' X3 4tt VI -z 2 v /l + y 

*,J = 2,3, — — = - — — [y 1 (x 2 )y 1 (x 3 )] 0)0 . (B.15) 
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Therefore, the integrals reduces to products of simple integrals of the kind 

h\e, m\ i, m, t, m"\ = J dx [Y e , m ,(x)]* Ye m (x)Y e „ m „(x) , (B.16) 

and 

I 2 [a, b, c, d, e, /, n, m] = J + * dx(l - x) a (l + x) b P^ d {x)P^{x) , (B.17) 

where the non-negative numerical coefficients a,. . .,/ are integers or half-integers and 
P£' d {x), P^i*{x) are Jacobi polynomials. Such integrals can be calculated analytically, 
or by using simple quadrature formulas (for more details, see the appendix of Ref. [85]). 

Eq. (B.6) can be used to evaluate the transformation coefficients of the HH functions 
with grand angular momentum G + 2 and n 3 > 0, once those corresponding to G are 
known. For the HH function with n 3 = one needs another recurrence formula, obtained 
from that one applied to the HH function with n 3 — 0, n 2 + 1. By proceeding as in the 
previous case, one obtains 

(p),G+2,L , . (p),G+2,L 

+ V a (p) ' G ' L 

^ a h/2,hM,n2fi-/'{/ 2 / 3 ,L' 2 \n' 2 \n' 3 ' 



X 



<Ag> + 2(6^-0^ , (B.18) 



where 



A _ • /V n 2 +1,0 O CR1Q\ 

" K-i 3 (l + «3 + / 9 3 )(2 + a 3 + / 93) ' 1 } 

n2 Aftf* 2(n 2 + \){n 2 + a 2 + fa + l)(2n 2 + a 2 + (3 2 ) ' 1 ' ) 

h , _ ^+i',o (2n 2 + a 2 + + l)(2n 2 + a 2 + & + 2) 

« 2 ^.JA 2 (n 2 + l)(n 2 + a 2 + / 3 2 + l) ' 1 ' } 



Kl'+ifi (n2 + a 2 )(n 2 + (3 2 )(2n 2 + a 2 + + 2) 



J^-ifl (n 2 + l)(n 2 + a 2 + (3 2 + l)(2n 2 + a 2 + (3 2 ) 

x (1 - 5 n2fi ) , (B.22) 

with a 2 = h + 1/2 and f3 2 = i 2 + 1/2. 

In practice, the recurrence relations are used as follows. After the calculation of 
the TC for the function with (n 2 = 0,n 3 = 0), those of the state (0, 1) are obtained by 
Eq. (B.6) and those of the state (1, 0) by Eq. (B.18). It can be noticed that, in the latter 
case, the term proportional to A n , 2 does not contribute, as c' n2=0 = 0. Then, Eq. (B.6) is 
used again for calculating the TC of the states (0, 2) and (1,1). The coefficients of the 
state (2, 0) can be now derived from Eq. (B.18), since those of the state (0, 2), entering 
that expression, are already available. The procedure can then be continued for larger 
values of (n 2 , n 3 ). 
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Appendix C. The Correlation Factors 

The correlation factors in Eqs. (83) and (85) have the Jastrow form, namely they are 
products of one-dimensional functions. Therefore, for each channel a few functions 
need to be chosen. One possibility would be to determine them by some preliminary 
variational procedure, as in Ref. [130] for a simplified problem (central interaction). 
However, such an approach would be numerically very involved and the following simpler 
procedure is preferred. 

In a generic nuclear system, when all the remaining particles are far from a given 
pair, the dependence of the total wave function on the coordinates of the pair of particles 
is mainly determined by their mutual interaction. Therefore, the radial wave function 
of the relative motion of the i, j pair in the angular-spin state (5 = 2Sl3+1 (£p)j , can be 
approximately described by the solution of an equation of the form 



where Vppi(r) =< /3'\V(i,j)\(3 > and V(i,j) is the interparticle potential. Depending 
on the quantum numbers, the state f3 can be a single state or coupled to other ones. 
The additional term Xpp'(r) in Eq. (C.l) has the role of simulating the average effect 
on the pair from the other particles. There is a large arbitrariness in choosing Xpp'(r), 
since the important condition to be satisfied is |A/3/3/(r)| <C |Vfo(?/(r)| when r is small. In 
Ref. [92] , it was found that a satisfactory choice is 



The value of 1/c should be greater than the range of the potential Vpp>(r), but its precise 
value has been found to be unimportant (the value adopted is 1/c = 2.0 fm). The depth 
Ap is then fixed so that (f>p(r) satisfy some appropriate healing condition. In this work, 
we have chosen <f>p(r) — > r^ 3 when r — > oo. The functions /3 (r)/r^ have been used to 
construct the appropriate correlation factor F a for a given channel a. 

Let us consider first the A = 3 case. The functions f a (r) are related to the reference 
pair (i,j) characterized by definite values of the angular momentum, spin and isospin for 
each channel; therefore these functions can be taken as solutions of Eq. (C.l). However, 
since the total wave function has been constructed in the LS'-coupling scheme, in general 
the total angular momentum j p of the reference pair has not a definite value. For the 
first three channels of Table 1, such a problem does not exist; in fact, l\ — and, since 
J =1/2, one has j p — S aa . As a consequence, the correlations f a , ce — 1, 3, correspond 
to the states 3 Si, 1 S and 3 Di, respectively. The channels with a > 3 have a minor 
relevance to produce the structure of the system than the first three ones. Therefore, 
we have chosen f a {r) = 03 5l (r) for the channels with £ 2 = 0, S aa = 1, T aa = 0, 
fa(r) = 0i5 o ( r ) for the channels with £ 2 = 0, S aa = 0, T aa = 1, and f a (r) = <fis Dl (r)/r 2 
for the channels with £ 2 = 2, S aa = 1, T aa = 0. Furthermore, for the channels with 
£ 2 = 1, we have used the solution of Eq. (C.l) by taking into account only the central 
part of the pair potential in the state 3 P. Finally, for the PHH expansion, g a {r) = 1, 




h 2 r d 2 2d £p(£p + 1) 



6pp>+Vpp,(r)+\pp,(r)}<f>p,(r) = ,(C.l) 



m idr 2 r dr r 2 



A/3/3' (r) = Apexp(-cr)5pp> . 



(C.2) 
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Table CI. Correlation functions for the 22 channels used for the CHH 
expansion of the ground state of the four-nucleon system. The functions 
fx = <f>s Sl , f 2 = <f>is , /3 = 0s Dl /r 2 , and / 4 = fopjr are the solutions of 
Eq. (C.l) for the state (3 = 2S+1 Lj. The function <p is calculated as explained 
in the text. 
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while in the CHH case, the functions g a (r) have be taken equal to 0(r), the latter 
function being chosen to be the solution of Eq. (C.l) with £p = and 

VpAr) = \(Vs C) (r) + V^ c) (r)) V • (C.3) 

Here Vg°\r) (V^ c \r)) is the central part of the nuclear potential projected on the state 
1 S'o ( 3 5i)- At large interparticle distances, the correlation factors go to 1 in order to 
recover the HH expansion, which is well suited for describing such configurations. 

Let us now consider the A = 4 case. The four-body correlation factors are given 
in term of the functions f a , g a and h a in Eq. (85). The functions f a ( r ij) related to the 
reference pair have been determined by using the same criteria as for the A = 3 
case. The functions g a {r- ik ) and h a (r km ) with k, m different from the reference pair 
indices i, j (see Eq. (85)) correlate pairs which are not in a definite angular-spin-isospin 
state. For simplicity reasons, for the channels constructed with the set A of the Jacobi 
vectors (a = 1 — 8 and 12 — 19, as listed in Table 2), we have chosen g a = h a = <f), where 
the function has been calculated as in the A = 3 case, namely using the potential 
given in Eq. (C.3). Other choices for these functions have been proved not to influence 
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appreciably the final result. 

The correlation functions of the channels constructed with the set B of the Jacobi 
vectors (a = 9 — 11 and 20 — 22) have been chosen in a similar way. In this case, the 
functions correlating pairs with definite angular-spin-isospin quantum numbers are f a 
and h a , and they have been chosen to be (p3 Sl or <p3 Dl depending on the value of £3 and 
£1, respectively. The functions g a have been chosen to be 0, as in the other channels. A 
summary of the correlation factors used in the various channels is given in Table CI. 



Appendix D. The Matrix Elements of the Interaction in the HH Expansion 



In this appendix, a more detailed discussion how to compute the matrix elements of the 
interaction between the basis functions S^ Tj7r /;(p) of Eq. (77) is reported. The states 
S^ Tj7r are constructed in terms of the Jacobi vectors defined with permutation p — 1 
and the standard choice of Jacobi vectors given in Eq. (5). In this case, the Jacobi 
vectors are specified in Eq. (5), namely 



xjv = r 2 - ri , 



x ^v-i = A/ 3 



r 3 



ri + r 2 



(D.l) 



Let us remember the general expression of the jj-coupling spin-isospin-HH states S^ Tj7r : 



^GTJtt 



{ (Yi n (x n )S 2 ) _ (xjv-i)s 3 ) . ■■■[Y h {x 1 )s A ) } 



x <^ (T 2 t 3 ) T3 t 4 



t A 



(D.2) 



where S 2 (T 2 ) is the spin (isospin) state of the particles 1 and 2, etc. Moreover, let us 
use the short-cut {T} = {T 2 , T 3 , . . . , Ta-i, Ta}, where Tj is the total isospin of particles 
1, . . . ,i and Ta = T. Let us consider first the matrix elements of the NN interaction, 
which in general is assumed to be non-local and non-isospin conserving. In configuration 
space, the matrix element in Eq. (99) is given by 



v G'T'v',GTv 



< ~^' T ' J >(1,2)|£^> 



= / cf XatgPx^GTXa^i • • • cf Xi 



^G'T'Jtv, 



(tftf) fl'(p') 



(D.3) 



where the dependence of V on the spin and isospin operators of particles 1 and 
2 is understood. Above, Qn denotes the hyperangular variables constructed from 
{x 1; . . . , Xjv-i, Xjv} and Q' N those constructed from {x. 1: . . . , Xjv-i, Xjy}. Moreover 

p = \Jxi + \- x 2 N and p' = ^xi + h a;^ r _ 1 + {x' N ) 2 . Note that V cannot change 

the total angular momentum j N of particles 1, 2. Let us suppose to know 



V J 



j N ,{T'},{T} 



i' N S' 2 ,£ N S2 \ X Ni X N) 



dx' N alx N (y v {£ N )S'^\ 
x V(x' N ,x N ) 



JN 



T" 



T'T Z 
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x [Yt N (x N )Si 



3N 



(T2*3)t 3 *4 



T 4 



■t A 



(DA) 



TT, 



which is usually quite easy to obtain. 

The integration over g? 3 xjv_i • • -d 3 ^ is easily performed applying the properties of 
the HH function. First of all, the hyperspherical angles ip 2 , • • • , ty?jv-i do not depend on 
xn, and therefore are the same for both and Q' N . Second, it results that 

rf 3 XAT_i • • • d 3 xi = (p N - 1 ) D ~ 4: dp N - 1 d£L N - 1 , (D.5) 

where iljv-i = {^1, • • • , £jv-i, <^2, • • • ? Vjv-i} are the hyperspherical coordinates related 
to the Jacobi vectors {xi, . . . , Xjv-i} and 



Pn-i = H \-Xfi_-y ■ (D.6) 

Third, x N (x.' N ) depends only on p and (p N (p' and (p' N ) and hence V does not depend on 
the Qn-i variables. Therefore, we can integrate over o!Qn~i and use the orthonormality 
relation of Eq. (32). Also the traces over the spin variables of particles 3, . . . , A can be 
easily computed, and the final results is: 

J x 2 N alx N (x' N ) 2 dx' N (pN-i^^dpf 



yJ-rc 

v G'T'v',GTv 



>N-1 



'£' N S' 2 ,£ N S2 ( x Ni x n) 



where 



2,2 
X N ' Pn-1 ' 



(p' 



/\2 



[X 



n) 



+ Pn- 



1 5 



(D.7) 



(D.8) 



cosip N = x N /p, cos ip' N = x' N / p' , 

and Gn-i is given by Eq. (30). Finally, the three-dimensional integrations involved 
in Eq. (D.7) can be accurately performed with standard numerical techniques (Gauss 
integration) [82]. For local potentials, the functions V(x' N , xn) reduces to V(xn)5(x' n — 
x N )/x 2 N . 

For the 3N interaction W = W(x.' N , x' iV _ 1 , xjy, Xjv-i), the first step is to know the 
projection over three-body angular-spin states, namely 
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t'N S 23N e N-li'N-l> e N S MN e N-dN- 
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x 



x 



x 



(^2*3)7^4 
Ye N (x N )S 2 

(^2*3)t 3 *4 



rpl 



t 4 



t A 



Ye N (x' N )S' 2 ) f (Y^Jx'^ss) f 

t 

W(x5 v ,x5 v _ 1 ,x JV ,x JV _i) 

T'T Z 

y^ JV _ 1 (^_i)s 3 ) 

' Jiv-i J Jn-i 

■t A ) , 

J TT, 



(D.9) 
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where 



Then 

< ~s t,j *\w\~z tj * >= 

x 2 n cLxn x^ r _ 1 (ixjv_i (x' N ) 2 dx' N (x' N _ 1 ) 2 dx' N _ 1 {pn-2) D 7 dpN-2 

tttJ N—l ll^ 1 '})!^} (ft \ 

X fyiv-arfVa ' ' ' < *tt»ii < W-i,j; r _i ' ' ' ^2,J 2 

X ' ' ' ^1/5^-2^2 • • • <W 2 , (D.10) 

p 2 = X 2 N + X 2 N _ l + p 2 N _ 2 , COS(f N = X N /p, 

cos y?Ar_i = Xat_i/ x 2 N _ x + p^_ 2 , (D.ll) 
and analogously for the primed quantities. Moreover, 

N-2 

G N -2 = ^(4 + 2raO , 
i=i 

= 4_! + 2n'^_ 1 + GW_ 2 . (D.12) 
For momentum-space potentials, similar results are obtained, with the only change 
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